{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gl/venv0825/lib/python3.6/site-packages/ipykernel_launcher.py:46: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 26469 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 27874 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 26041 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 21521 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 24402 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 19968 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 21270 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 31354 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 38388 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 35889 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 26469 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 27874 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 26041 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 21521 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 24402 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 19968 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 21270 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 31354 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 38388 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 35889 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAk7UlEQVR4nO3de5hcdZ3n8fe3qvqW0LkgHSA3AxiY5SKgLTCuuCqo4CKIOooz6/3ZDPPA7Lq788gw7OMyOjzjZXT3cbyRecZdfZYRHRFBvADxMjq7C5ggAhHBQGBIJySBXLqT9K2qvvvHOae6ui7ddTvnVIfP63n66apzTp36plJ9vud3N3dHRESkXCbtAEREpPsoOYiISBUlBxERqaLkICIiVZQcRESkSi7tADrhuOOO83Xr1qUdhojIgrJly5bn3X2o1r6jIjmsW7eOzZs3px2GiMiCYmbP1NunaiUREami5CAiIlWUHEREpIqSg4iIVFFyEBGRKl2bHMzsEjN73My2mdmfpx2PiMiLSVcmBzPLAl8ELgVOB95jZqenG5WIyItHVyYH4Dxgm7s/5e5TwK3AFSnHJJIKd+fbW3YwMV1IOxR5EenW5LAKeLbs+Y5wW4mZbTCzzWa2ee/evYkGJ5KkbXsO8Wf/+Gt+9ri+55Kcbk0O83L3je4+7O7DQ0M1R3+LHBUm80UApgrFlCORF5NuTQ4jwJqy56vDbSIvOsVwtcZiUas2SnK6NTn8ElhvZieZWS9wFXBnyjGJpKIQJoWCkoMkqCsn3nP3vJldC9wNZIGvuvvWlMMSSUVUcihovXdJUFcmBwB3/wHwg7TjEElbvqCSgySvW6uVRCQUlRiUHCRJSg4iXa4YdlIqqlpJEqTkINLlVHKQNCg5iHS5onorSQqUHES6XD5MCqpWkiQpOYh0uajEkFfJQRKk5CDS5TRCWtKg5CDS5WZGSKcciLyoKDmIdDmNkJY0KDmIdLmo5KBqJUmSkoNIl4saolVykCQpOYh0OY1zkDQoOYh0OY2QljQoOYh0OZUcJA1KDiIxKhad//l/tnNkKt/yOQoaIS0pUHIQidETe8b4y+/9hp8/sbflc4TLOajkIInquuRgZjea2YiZPRT+vCXtmERaNZ0PLuhThdYv7IVwzm6VHCRJ3boS3H93979JOwiRdhU6MPVFNDI630aCEWlW15UcRI4m0V1/O5PmaYS0pKFbk8O1ZvawmX3VzJbXOsDMNpjZZjPbvHdv6/W5InGaWf+59YmRNEJa0pBKcjCzTWb2aI2fK4AvA6cA5wC7gM/WOoe7b3T3YXcfHhoaSi54kSbMjFFo4xylEdKdiEikMam0Obj7xY0cZ2Z/B9wVczgisZmZUbX17KApuyUNXVetZGYnlj29Eng0rVhE2pXvwEI9eQ2CkxR0Y2+lT5vZOYADTwN/nGo0Im0oFNq/sBe1EpykoOuSg7u/N+0YRDqlE/MiaYS0pKHrqpVEjiadWP9ZE+9JGpQcRGLUifaCokoOkgIlB5EYRb2U2kkOapCWNCg5iMQoGt/QVslB1UqSAiUHkRh1YvoMNUhLGpQcRGKU78CFvTTxnkoOkiAlB5EYlcYotDH3hUZISxqUHERilO/A9BkzcyspOUhylBxEYtSRcQ6lBNORkEQaouQgEqPOtDmoWkmSp+QgEqNCB9ocClrsR1Kg5CASo0IHR0hrnIMkSclBJEb5DjQma24lSYOSg0iMOjkITslBkqTkIBKj0vQZ7bQ5aIS0pEDJQSRGKjnIQpVKcjCzPzCzrWZWNLPhin3Xm9k2M3vczN6cRnwindKJrqylEdIqOUiC0loJ7lHg7cDN5RvN7HTgKuAMYCWwycxOdfdC8iGKtK8TS3x2YiCdSLNSKTm4+2Pu/niNXVcAt7r7pLtvB7YB5yUbnUjndGT6jDAnqFpJktRtbQ6rgGfLnu8It1Uxsw1mttnMNu/duzeR4ESa1Yn2giixaIS0JCm2aiUz2wScUGPXDe5+R7vnd/eNwEaA4eFh/dVIV+rEKm6lHk9qc5AExZYc3P3iFl42Aqwpe7463CayIHWizaG0hrQm3pMEdVu10p3AVWbWZ2YnAeuBB1KOSaRlHSk5aG4lSUFaXVmvNLMdwO8D3zezuwHcfSvwLeA3wI+Aa9RTSRayTs+t5EoQkpBUurK6++3A7XX23QTclGxEIvHIh3VBnSg5ABQdstZ2WCLz6rZqJZGjSifWfy6f7lvdWSUpSg4iMepEN9TirJKDkoMkQ8lBJEb5Do6QrnwsEiclB5EYFTuwFkN5aUE9liQpSg4iMYraC/LtTJ9RXnJoY+pvkWYoOYjEaKYra3vnyGaCLkoqOUhSlBxEYtSRifeKTm82+FPV/EqSFCUHkRhF7QVtNUi705vLlB6LJEHJQSRGUZtDW11Zi9ATlhzUW0mSouQgEqNOLNRTcKc3HBatyfckKUoOIjFqd/oMd6dQdHrCaqV2ej2JNEPJQSRGUU5ota0gen1UraQR0pIUJQeRGEV3+u6ttTtEJY7eUptD52ITmYuSg0iMygettdLuEJUUSr2V1CAtCVFyEIlReXVSKxf2ypKDqpUkKU0nBzNbbGbZdt7UzP7AzLaaWdHMhsu2rzOzcTN7KPz5SjvvI5K2WVNftHBhj17Tk7Oq84nEad7FfswsA1wF/BHwKmAS6DOz54HvAze7+7Ym3/dR4O3AzTX2Penu5zR5PpGulG9zXqSonSJqkG6nS6xIMxopOfwUOAW4HjjB3de4+wrgNcB9wKfM7N8186bu/pi7P950tCILTKHg5MJ5kVrphpqvSA6qVpKkNLJM6MXuPl250d33AbcBt5lZTwdjOsnMfgWMAv/V3X/RwXOLJKrgTl8uQ36q0FK1UlRyUIO0JG3e5FCeGMxsObASGAeedvdi5TFlx24CTqhxyhvc/Y46b7cLWOvuL5jZK4HvmtkZ7j5a4/wbgA0Aa9eune+fIZKKfNFZ3Jvl8FShtQbpMKH0aeI9SVgjbQ5LgWuA9wC9wF6gHzjezO4DvuTuP618nbtf3Gww7j5J0KaBu28xsyeBU4HNNY7dCGwEGB4e1l+MdKVCcWbSvHwLbQ6FimolTbwnSWmkWunbwNeBC939QPmO8O7+vWZ2srv/fbvBmNkQsM/dC2Z2MrAeeKrd84qkIZr6oi8XdO5rpeQQNVOot5IkrZFqpTfOsW8LsKXZNzWzK4G/BYaA75vZQ+7+ZuC1wMfNbBooAleHbRsiC050He9rY7rtUldWzcoqCWuk5ICZ5YCCu7uZrQHOJ+hy+qtW3tTdbwdur7H9NoJGbpEFL+qd1E5jcrRIUK+SgyRs3q6sZvbvgT3AM+HjHwPvBG41s+tijk9kwSpU9DRqrc2BWedQV1ZJSiMlh48QjHMYBB4DXuruz5vZIuCXwKfiC09k4YqSQ18bF3ZNvCdpaSQ5TLn7fmC/mW1z9+cB3P2ImU3FG57IwjVTcggapNuZeK9Hy4RKwhpJDgNmdi5BFVRv+NjCn/44gxNZyPJVd/3N3/ZXdmXVOAdJSiPJ4TngczUeR89FpIZStVJPG20O0ZTd2WgKDiUHSUYjXVlfl0AcIkedyjaHlrqyquQgKWlkhPTb59rv7t/pXDgiR4+q5NDOeg5qc5CENVKt9Nbw9wrg1cBPwuevB/4voOQgUkNlm0NLDdKaeE9S0ki10gcBzOwe4HR33xU+PxH4X7FGJ7KARQ3QfT1Bb6WW1pCuGCGtcQ6SlGZWglsTJYbQbkDToYrUURrA1kbJoXqcg5KDJKOh6TNCPzazu4FvhM/fDWzqfEgiR4fOTJ+huZUkHQ0nB3e/Npww77Xhpo3hHEkiUkMnG6R7spqVVZLVSG8lcw8qOutNmFd+jIgEOpEcojYG9VaSpDW0hrSZ/amZzWpfMLNeM3uDmX0NeH884YksXJ2YPqNq4j2VHCQhjVQrXQJ8CPhGuADPfmCAILHcA/yPVqfuFjma5au6obYwfYZXNkh3KDiReTTSlXUC+BLwJTPrAY4DxitXhROR2aqrlZo/R1RSyGbCNgdVK0lCGl3s52M1tpU/3ePuX2n0Tc3sMwSD66aAJ4EPRsnGzK4HPgwUgP/g7nc3el6RblI1urmFkkNU+shlMmQzpmolSUyjvZUuAK4imIm1lq8BDScH4F7genfPm9mngOuB68zs9PB9zgBWApvM7FR3LzRxbpGuUFlyaGeEdCYDWTNNvCeJaXQQXMHdR939YK0foKlvrLvf4+758Ol9wOrw8RXAre4+6e7bgW3Aec2cW6RbVLc5tD5COpuxoOSgaiVJSKPJYb5vZDvf2A8BPwwfrwKeLdu3I9xWxcw2mNlmM9u8d+/eNt5eJB6l6TPC3krtjHPIWpAcNM5BktJotVKPmS2ps8+AbNVGs03ACTWOv8Hd7wiPuQHIA7c0GEeJu28ENgIMDw/rL0a6TtQA3Va1kkfVSkbGNAhOktNocriPYC3pen5YucHdL57rhGb2AeAy4KKyAXQjwJqyw1aH20QWnE5OnxGVHFStJElpZuI9m+OnKWZ2CfBR4HJ3P1K2607gKjPrM7OTgPXAA82eX6QbdGJepFJyyKpaSZLVaMnhfDrbW+kLQB9wb9gl9j53v9rdt5rZt4DfEFQ3XaOeSrJQzXRDbf3CXl5yyJiSgySn0eRQcPfRejvNrNneSi+bY99NwE3NnE+kG0XdUHPhXX9L02eU9VbKqeQgCeqG3koiR6V82V1/rsX2gtI4BzMyGdMIaUlMbL2VRF7sCmVTX2TNyBdan3ivNM5BJQdJSGy9lURe7ArlU19krbWJ98LXZCwogbSQX0Ra0sxKcE33ShJ5MSvvaZRro80hmzEsrFZSyUGSklZvJZGjXr4DYxQKxeD10XnyLZQ+RFqRSm8lkReDqEqonTaHojuZsNtI0B22kxGK1KfeSiIxiS7kuYyFbQ6tjXMolRw0QloSpN5KIjEpNSZnjFwm01I31ELRyYQL/WQ0zkESpN5KIjHJF51cdGG3VteQnjlH1lDJQRKj3koiMSkUvbS8Zy6TodDKOAefOYfmVpIkqbeSSExm3fW32JW1WHQyFpU+tBKcJEe9lURiki9rL8hlW+3KWlb6yBqT0+quJMlQbyWRmBRmtTm0PgiuvOSguZUkKeqtJBKToL0guP8KZlRt/q6/WJzd5qAR0pIU9VYSiUmhMLvNoZXG5Pys3koqOUhy1FtJJCb5irv+6RaGNwcjpMvHOXQ0RJG6UumtZGafAd4KTAFPAh909wNmtg54DHg8PPQ+d7+60fOKdJNCsTgrOYxPtzlC2lqrmhJpRVq9le4Frnf3vJl9CrgeuC7c96S7n9Pk+US6TsEpVQnlWmwvKBQplRxanYJDpBWp9FZy93vcPR8+vQ9Y3czrRRaC2SWHTGvjHNzJRhPvmaHcIElpNDn0mNmSOj9Laa+30oeY3aB9kpn9ysz+ycwurPciM9tgZpvNbPPevXvbeHuReOQL5W0OtNwgXT7xnkoOkpRO9FYyavRWMrNNwAk1jr/B3e8Ij7kByAO3hPt2AWvd/QUzeyXwXTM7o1aVlrtvBDYCDA8P6y9Guk7l9BmtjpDOlo2VUHKQpMTWIO3uF891QjP7AHAZcJF70D/P3SeByfDxFjN7EjgV2NxgnCJdo+Czu7K21uYwu/ShifckKak0SJvZJcBHgX/j7kfKtg8B+9y9YGYnA+uBp5o5t0i3KJ9uu9W5lcpHSLfabiHSikaTQ6enz/gC0Afca8EXP+qy+lrg42Y2DRSBq919X5PnFukK+Q4MgisWnd5cJjwHGiEtiUll+gx3f1md7bcBtzVzLpFuVT7ddq7F5DBrym6NkJYENdsgXa/N4UcdiUbkKFIoOr09wX1Ty9VK5Q3S6q0kCWooObj7X8YdiMjRJl/R5tDK6ObKEdKqVpKkNDrOQUSaVCgW217sp7JRW9VKkhQlB5GYFIrManNo5a6/6BoEJ+lQchCJSXnJIdOBNgclB0mSkoNITGYtE9pqV1afmXgvE86t5KpakgQoOYjEpHyZ0Gwm01J7Qb6i3QLQ5HuSCCUHkZiU9zTKZQz35gexTeWL9GRnJ4e81nSQBCg5iMRkfKrAQG8wziEXXuCnmlzKbWwiz2B/DwCLwnMdmSx0MEqR2pQcRGLg7oxOTJcu7IN9wZCiQ5P5uV42S75Q5MhUgSXhOaLfYxONn0OkVUoOIjGYzBeZLjhLBoKksGQguLCPjk83fI4oCUTnGOwPfo9ONH4OkVYpOYjEIEoCgxV3/aNN3PVHSaB0jijBKDlIApQcRGIQJYEl/VHJIbzrb6Xk0F9RchhXtZLET8lBJAZj4d19ZXtBM3f9USKJSgytnEOkVUoOIjGISg6D/ZVtDq1UK80+hxqkJQlKDiIxKJUc2rjrn6maCl57TF/zVVMirUotOZjZJ8zsYTN7yMzuMbOV4XYzs8+b2bZw/yvSilGkVVEJIbrr7+/J0JO1pi7sldVK2Ywx2JdTyUESkWbJ4TPu/nJ3Pwe4C/hYuP1SgrWj1wMbgC+nE55I68YqehqZGUv6e1oqOUQlhuB8ObU5SCJSSw7uPlr2dDEz61BfAXzdA/cBy8zsxMQDFGnD2ESejMHi3pkVdJcM9DTV5jA2Mc1gX640bUZ0jjElB0lAo8uExsLMbgLeBxwEXh9uXgU8W3bYjnDbrorXbiAoWbB27drYYxVpRjQ62qzswt7kXf/oeL5UpRQZ7M+pK6skItaSg5ltMrNHa/xcAeDuN7j7GuAW4Npmzu3uG9192N2Hh4aG4ghfpGVjE/nS2IZIUHJoplpputRmUTpHfw9jkyo5SPxiLTm4+8UNHnoL8APgvwEjwJqyfavDbSILxuj4NIN91Xf9uw5ONHyOsYnpUk+l8nP8bo9KDhK/NHsrrS97egXw2/DxncD7wl5LFwAH3X1X1QlEulgwm2r1XX9zvZXqlD7U5iAJSLPN4ZNmdhpQBJ4Brg63/wB4C7ANOAJ8MJ3wRFo3OjHNmmMXzdrW7IV9dGKa0/oHZ20b7A+6srr7rPYMkU5LLTm4+zvqbHfgmoTDEemo2iWHHBPTRabyRXpz8xfaxybypXmVZs7RQ6HoHJkqsLgv1f4kcpTTCGmRGIzWaC+Ymf5i/tJDsehBm0NVbyVNoSHJUHIQ6bBi0Tk0WfuuHxqbtvvwVJ6iU136GNCaDpIMJQeRDjs0lcd95i4/0sy03WMV8ypFZkoOSg4SLyUHkQ6rXMEt0szke6MVE/fNnENrOkgylBxEOqxyFbhIM9N2V07cFxnUmg6SECUHkQ6rVyXUTMmhcrGg0jlKbQ4qOUi8lBxEOmym5FCnMbmBNof61UpR6UMlB4mXkoNIh0VzH1Umh4GeLLmMNdbmUKdaqS+XoTebUVdWiZ2Sg0iHzTRIz77rN7OGp+0em6idYMxMazpIIpQcRDqsXrUSND5t9+hEnr5chr5ctmpfsKaDSg4SLyUHkQ4bm+fCfrCB9oK9Y5Mcd0xfzX1L+nMNnUOkHUoOIh32wuEpli/qrblvxWA/zzUwbffI/nFWLRuouW9osJ89o41P/S3SCiUHkQ7bdXCcE5f119y3evkAI/vHCeaXrG/kwDirltdODquW9TNyYLztOEXmouQg0mE7D0ywss5d/+rlA4xN5udslJ4uFNl1cJzVdZLDicsGGJvIawoNiZWSg0gHuTs7D4yzcmn9kgPAs/uP1D3HcwcnKDp1q5WixNPMqnIizUolOZjZJ8zsYTN7yMzuMbOV4fbXmdnBcPtDZvaxNOITadW+w1NM5otzlByCBYB27K9fLRRVGc1VrVR+nEgc0io5fMbdX+7u5wB3AeVJ4Bfufk748/F0whNpzc4Dwd38XNVKADvmKDlEiSNKJJWic+9UcpAYpZIc3H207OliYO7WOZEFYufB4IK9cmnt5LB0oIdj+nJzlxzCfSfWqZpaMdhPNmNKDhKr1NYZNLObgPcBB4HXl+36fTP7NbAT+DN335pGfCKtiC7YK+v0VjIzVi8fmKda6QhDg33091SPkwDIZowTlvSXSikicYit5GBmm8zs0Ro/VwC4+w3uvga4Bbg2fNmDwEvd/Wzgb4HvznH+DWa22cw27927N65/hkhTdh4Ypy+X4djFtcc5AGFymLtaqV5PpcjKZf0qOUisYksO7n6xu59Z4+eOikNvAd4RvmbU3Q+Fj38A9JjZcXXOv9Hdh919eGhoKK5/hkhTdh4MurGaWd1jVi9fNOdYh5ED9QfARVYuGyhVYYnEIa3eSuvLnl4B/DbcfoKFf1Vmdh5BfC8kH6FIa3YeGK9bpRSZa6xDsRh0ha3XUymyctkAzx2coFBUc53EI602h0+a2WlAEXgGuDrc/k7gT8wsD4wDV/l8Q0lFusiuAxNcuL5mYbekfKzD0kVLZ+3bMzbJdMHr9lSKrFw2wHTBef7QJMcvmTsZibQileTg7u+os/0LwBcSDkekI6YLRXaPTXDiPFVC5WMdzlw1OzmMHAjaIlbPV620dGasg5KDxEEjpEU6ZPfoBO4zg9TqWXNskBye3Huoat/DOw4CcOoJg3OeozRKWj2WJCZKDiIdMjM+Ye67/qUDPaxfcQwPbN9Xte+XT+9j1bKBeRuko6qpp1843GK0InNTchDpkMd3jwHwshXHzHvsBSe/hM1P7yNfKJa2uTsPbN/Pq9Ytn/f1g/09rD12EVt3Hmw9YJE5KDmIdMgjOw7yksW9dUc2lzv/5GM5PFXg0Z0zkwU8/cIRnj80yatOOrah9ztr1VIeGVFykHgoOYh0yCMjBzlz1dI5xzhEzgsTwP1PzfTUfmB78Pj8BpPDGauW8Oy+cQ4cmWohWpG5KTmIdMDEdIHf7TnEWRW9j+pZMdjPKUOLub+s3eGB7fs5dnEvpwzNXy0FlN5r687ReY4UaZ6Sg0gHPLZrlELROXPVkoZfc/7JL+GX24N2h0LRuX/7Cwy/dHlDJQ+AM1cGyUFVSxIHJQeRDng0vEBXjluYyxtOW8HYZJ7P3fsEX/zpNnbsH+eys1c2/Prli3tZtWyg9N4inZTarKwiR5NHRg6yfFHPvF1Qy130r1bwnvPW8qWfPYkZvO2clbz15Sc29b5nrVqq5CCxUMlBpAMeGRltuDE6Ymb85eVncOH64zjt+EH+6sqzmno9wFmrl/L0C0c4eETrSUtnqeQg0qb9h6d4YvcYF/3eKU2/tjeX4esfOo9C0cllm79Xu+DkoGfTTx7fzZXnrm769SL1qOQg0qYfPvochaJzyZkntPR6M2spMQCcu2Y5K5f2c9evd7X0epF6lBxE2nTnr0c4eWgxZ6xsvKdSp2QyxmVnr+Tnv9urqiXpKCUHkTY8d3CC+7fv4/KzVzbdXtApl738RKYLzt1bn0vl/eXopOQg0oa7Ht6JO1zeRBfUTjtr1VLWHruI2x7cUXd1OZFmKTmItGjf4Slu/vlTnLt2GSc3OKo5DmbGB169jvu37+OOh3amFoccXVJPDmb2X8zMo7WiLfB5M9tmZg+b2SvSjlGkkrtz/Xce5uCRaW5621lph8P7X72OV6xdxo3f28qeMa3xIO1LtSurma0B3gT8S9nmS4H14c/5wJfD3yJNmy4UmcwXmZguMDFdYCpfJJfJ0NeToS+XoS+XpTeXYWximudGJ9i25xC/232IbXsO8cLhSSamixx3TB8vfckizly1hDXLFzGVL/K1//c0d2/dzfWX/h6np9AQXSmbMT79zrN5y+d/wb/9/D/zkYvXc8rQMRyZyvPE7kM8sXuM7c8fZrpQpCebYeXSAVYtHyitHbFqefCzpL8Hd6foUCg6RXcKRafgTrHouENfT4b+XJZMJp02FkmGpVlHaWbfBj4B3AEMu/vzZnYz8DN3/0Z4zOPA69y9bl+94eFh37x5c9Pv/9iuUa7+31tq7qv3sTi1d9Q9voWPt97/Sb1TNRvr3K9p7vj6r2jlPer/u4vFmQtWdLFyd7IZw8zIGGTNyJhhFhw3kQ/mLGpWxmDtsYtYsaSfvlyGPaOTPLPvMBPTM2svDPbl+PCFJ/Gnb1hPtosukr9+9gA3fm8rv/qXA7O2H7+kj1OGjmGgJ8tEvsDOAxOMHBhnKl+cdVzGoNGPrC+XYaA3y0BPlkzYGB/9H5afovy/NfpOzt5W69h6Qcx81lH7f/mnP7NtvuO65/+sXGVY5c8Nq7nv9aet4MbLz2jx/WyLuw/X2pdaycHMrgBG3P3XFf9Rq4Bny57vCLfNSg5mtgHYALB27dqWYljUm+XcNcvmirH29rovqLe5/hex3ne03ivqH9/CezT5D2k21rlf09x7ZMILfy4b/M5mgnMUw7vcYnhnGz3OZax0h9vfky097s1lKBSdyXyRyXwh+D1d5Jj+HCsGgwvoyUOL6e/Jznr/fKHIk3sPs3t0gulCkeF1x7J0oGeOf3k6zl6zjO/8yat5ZOQghybz9OUyvGxokKWLqmMtFp3nD08ysn+ckQPjjOwfZ2wiTyZjZMPPeObxTOKdyhcZny4EP1PBT/mlPPovrHdhK12sZ/1fz30xh/mTSLSt0WTkPtffQDJqxVq2odbD8HUzW046bnHnAyPmkoOZbQJqjQy6AfgL4E3uftDMnmam5HAX8El3/+fwHD8GrnP3ukWDVksOIiIvZqmVHNz94joBnQWcBESlhtXAg2Z2HjACrCk7fHW4TUREEpJKbyV3f8TdV7j7OndfR1B19Ap3fw64E3hf2GvpAuDgXO0NIiLSed048d4PgLcA24AjwAfTDUdE5MWnK5JDWHqIHjtwTXrRiIhI6oPgRESk+yg5iIhIFSUHERGpouQgIiJVUp0+o1PMbC/wTNpxzOM44Pm0g2jQQol1ocQJCyfWhRInLJxYuznOl7r7UK0dR0VyWAjMbHO9kYjdZqHEulDihIUT60KJExZOrAslzkqqVhIRkSpKDiIiUkXJITkb0w6gCQsl1oUSJyycWBdKnLBwYl0occ6iNgcREamikoOIiFRRchARkSpKDjEzs2+a2UPhz9Nm9lC4fZ2ZjZft+0rKcd5oZiNl8bylbN/1ZrbNzB43szenGWcYz2fM7Ldm9rCZ3W5my8LtXfWZhjFdEn5u28zsz9OOJ2Jma8zsp2b2GzPbamb/Mdxe93uQpvBv55Ewps3htmPN7F4z+134e3kXxHla2Wf3kJmNmtlHuvVznYvaHBJkZp8lWJ/i42a2DrjL3c9MOSwguCgAh9z9byq2nw58AzgPWAlsAk5190LiQc7E9CbgJ+6eN7NPAbj7dV34mWaBJ4A3EqxZ8kvgPe7+m1QDA8zsROBEd3/QzAaBLcDbgHdR43uQtvLVIsu2fRrY5+6fDBPvcne/Lq0YK4X//yPA+QRLD3Td5zoXlRwSYsGSd+8iuNAuJFcAt7r7pLtvJ1hn47w0A3L3e9w9Hz69j2C1wG50HrDN3Z9y9yngVoLPM3XuvsvdHwwfjwGPEazVvpBcAXwtfPw1guTWTS4CnnT3bp+9oSYlh+RcCOx299+VbTvJzH5lZv9kZhemFViZa8Oqmq+WFdFXAc+WHbOD7rqIfAj4YdnzbvpMu/2zA4LqOOBc4P5wU63vQdocuMfMtpjZhnDb8WWrRD4HHJ9OaHVdxeybwW78XOtScugAM9tkZo/W+Cm/S3wPs78ou4C17n4u8J+BfzCzJSnG+WXgFOCcMLbPxhlLm7FGx9wA5IFbwk2Jf6YLnZkdA9wGfMTdR+my70GZ17j7K4BLgWvM7LXlO8NFwrqmjtzMeoHLgX8MN3Xr51pXV6wEt9C5+8Vz7TezHPB24JVlr5kEJsPHW8zsSeBUYHNacUbM7O+Au8KnI8Cast2rw22xauAz/QBwGXBReGFI5TOdRyqfXaPMrIcgMdzi7t8BcPfdZfvLvwepcveR8PceM7udoMput5md6O67wjaUPakGOdulwIPR59mtn+tcVHJIxsXAb919R7TBzIbCBivM7GRgPfBUSvFFDZSRK4FHw8d3AleZWZ+ZnUQQ5wNJx1fOzC4BPgpc7u5HyrZ31WdK0AC93sxOCu8kryL4PFMXtoH9PfCYu3+ubHu970FqzGxx2GiOmS0G3kQQ153A+8PD3g/ckU6ENc2qKejGz3U+Kjkko7LuEeC1wMfNbBooAle7+77EI5vxaTM7h6Bo/jTwxwDuvtXMvgX8hqAK55o0eyqFvgD0AfcG1zjuc/er6bLPNOxNdS1wN5AFvuruW9OKp8K/Bt4LPGJh92rgL4D31PoepOx44Pbw/zoH/IO7/8jMfgl8y8w+TDBl/7tSjLEkTGBvZPZnV/Pvq5upK6uIiFRRtZKIiFRRchARkSpKDiIiUkXJQUREqig5iIhIFSUHkZSY2UA4zUc0NuM/mdmDZvbu8Hmvmf08HEQpkih96UTqCGeqvYBgfAcEfy/31dlGM9vd/UaCeaG+4+6FcBqLVxGM/L0N+Ka7T5nZj4F3MzNFiEgilBxE5naVux8AsGDdiI/U2Vbv2Lm2/xHwh+FjC39XDjz6LvDXKDlIwlStJJKCcDqNk939aShNm/0IwTxQ3yw79FGCEoVIolRyEEnHccCB8g3u/tcEpYTybQUzmzKzwTCBiCRCJQeRdIwD/Q0e2wdMxBiLSBUlB5EUuPt+IGtmcyYIM3sJ8Ly7TycTmUhAyUEkPfcAr5nnmNcD308gFpFZlBxE0vNFZtYjqOcPgZsTiEVkFjVIi9S3B/i6mRXD5xngR3W20ex2d3/QzH5qZtlaa2SEPZq+6+5PdPRfJdIArecgIiJVVK0kIiJVlBxERKSKkoOIiFRRchARkSpKDiIiUuX/Axh02fBHBFnvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "def generate_array_data(M, N, d, wavelength, SNR, DOA, add_noise=1):\n",
    "\n",
    "    K = len(DOA.reshape(-1))\n",
    "    array_geom = np.expand_dims(np.array(np.arange(M)), axis=-1) * d\n",
    "\n",
    "    signal_i = 10 ** (SNR / 20) * (np.random.randn(K, N) + 1j * np.random.randn(K, N))\n",
    "    noise = np.random.randn(M, N) + 1j * np.random.randn(M, N)\n",
    "\n",
    "    phase_shift = 2 * np.pi * array_geom / wavelength * np.sin(DOA / 180 * np.pi)\n",
    "    a_i = np.exp(-1j*phase_shift)\n",
    "    array_signal = np.matmul(a_i, signal_i)  # 包络*相移\n",
    "    if add_noise:\n",
    "        array_data = array_signal + noise\n",
    "    else:\n",
    "        array_data = array_signal\n",
    "    array_covariance = 1 / N * (np.matmul(array_data, np.matrix.getH(array_data)))\n",
    "    return array_data,array_covariance\n",
    "def generate_steering_vector(M, d, wavelength, DOA):\n",
    "    K = len(DOA.reshape(-1))\n",
    "    array_geom = np.expand_dims(np.array(np.arange(M)), axis=-1) * d\n",
    "    phase_shift = 2 * np.pi * array_geom / wavelength * np.sin(DOA / 180 * np.pi)\n",
    "    a_i = np.exp(-1j*phase_shift)\n",
    "    return a_i\n",
    "# array signal parameters\n",
    "fc = 77e9     # carrier frequency\n",
    "c = 3e8      # light speed\n",
    "M = 16        # array sensor number\n",
    "N = 400       # snapshot number\n",
    "wavelength = c / fc  # signal wavelength\n",
    "d = 0.5 * wavelength  # inter-sensor distance\n",
    "SNR = 10\n",
    "DOA = np.array([-30,10])\n",
    "\n",
    "array_data,array_covariance = generate_array_data(M, N, d, wavelength, SNR, DOA)\n",
    "K = len(DOA)\n",
    "w,V = np.linalg.eig(array_covariance)\n",
    "w_index_order = np.argsort(w)\n",
    "V_noise = V[:,w_index_order[0:-K]]\n",
    "noise_subspace = (np.matmul(V_noise, np.matrix.getH(V_noise)))\n",
    "doa_search = np.arange(-90,90)\n",
    "p_music = np.zeros(len(doa_search))\n",
    "for doa_index in range(len(doa_search)):\n",
    "    a = generate_steering_vector(M, d, wavelength, doa_search[doa_index])\n",
    "    p_music[doa_index]=1/np.matmul(np.matmul(a.T.conjugate(),noise_subspace),a).reshape(-1)[0]\n",
    "p_music = p_music/np.max(p_music)\n",
    "p_music = 10*np.log10(p_music)\n",
    "plt.figure()\n",
    "plt.plot(doa_search,p_music)\n",
    "plt.ylabel(\"归一化空间谱(dB))\", fontproperties=\"SimSun\")\n",
    "plt.xlabel(\"来波方向(°)\", fontproperties=\"SimSun\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " def demo(self):\n",
    "        print(scipy.fft.fft(np.exp(2j * np.pi * np.arange(8) / 8)))\n",
    "        # Create two random square matrices A and B\n",
    "        A = np.random.rand(3, 3)\n",
    "        B = np.random.rand(3, 3)\n",
    "        # Create a vector x\n",
    "        x = np.random.rand(3)\n",
    "        # Demonstrate various operations\n",
    "        print(\"Matrix A:\\n\", A)\n",
    "        print(\"Matrix B:\\n\", B)\n",
    "        print(\"Vector x:\\n\", x)\n",
    "        # Addition\n",
    "        C = A + B\n",
    "        print(\"\\nAddition (A + B):\\n\", C)\n",
    "        # Subtraction\n",
    "        D = A - B\n",
    "        print(\"\\nSubtraction (A - B):\\n\", D)\n",
    "        # Matrix Multiplication\n",
    "        E = np.dot(A, B)\n",
    "        print(\"\\nMatrix Multiplication (A * B):\\n\", E)\n",
    "        # Element-wise Multiplication (Hadamard product)\n",
    "        F = np.multiply(A, B)\n",
    "        print(\"\\nElement-wise Multiplication (A .* B):\\n\", F)\n",
    "        # Inversion of matrix A\n",
    "        if np.linalg.det(A) != 0:  # Check if A is invertible\n",
    "            A_inv = inv(A)\n",
    "            print(\"\\nInversion of A:\\n\", A_inv)\n",
    "        else:\n",
    "            print(\"\\nMatrix A is not invertible.\")\n",
    "        # Eigenvalues and eigenvectors of matrix A\n",
    "        eigenvalues, eigenvectors = eig(A)\n",
    "        print(\"\\nEigenvalues of A:\\n\", eigenvalues)\n",
    "        print(\"\\nEigenvectors of A:\\n\", eigenvectors)\n",
    "        # Singular Value Decomposition of matrix A\n",
    "        U, s, Vt = svd(A, full_matrices=False)\n",
    "        print(\"\\nSingular values of A:\\n\", s)\n",
    "        print(\"\\nLeft singular vectors of A:\\n\", U)\n",
    "        print(\"\\nRight singular vectors of A:\\n\", Vt)\n",
    "        # Selection\n",
    "        # Select the first row of matrix A\n",
    "        first_row = A[0, :]\n",
    "        print(\"\\nFirst row of A:\\n\", first_row)\n",
    "        # Select the second column of matrix A\n",
    "        second_column = A[:, 1]\n",
    "        print(\"\\nSecond column of A:\\n\", second_column)\n",
    "        # Select a submatrix from A consisting of rows 1 and 2 and columns 1 and 2\n",
    "        submatrix = A[1:3, 1:3]\n",
    "        print(\"\\nSubmatrix of A (rows 1 and 2, columns 1 and 2):\\n\", submatrix)\n",
    "        \n",
    "        # Create two random complex matrices A and B\n",
    "        A = np.random.rand(3, 3) + 1j * np.random.rand(3, 3)\n",
    "        B = np.random.rand(3, 3) + 1j * np.random.rand(3, 3)\n",
    "        # Create a random complex vector x\n",
    "        x = np.random.rand(3) + 1j * np.random.rand(3)\n",
    "        # Demonstrate various operations\n",
    "        print(\"Matrix A:\\n\", A)\n",
    "        print(\"Matrix B:\\n\", B)\n",
    "        print(\"Vector x:\\n\", x)\n",
    "        # Addition\n",
    "        C = A + B\n",
    "        print(\"\\nAddition (A + B):\\n\", C)\n",
    "        # Subtraction\n",
    "        D = A - B\n",
    "        print(\"\\nSubtraction (A - B):\\n\", D)\n",
    "        # Matrix Multiplication\n",
    "        E = np.dot(A, B)\n",
    "        print(\"\\nMatrix Multiplication (A * B):\\n\", E)\n",
    "        # Element-wise Multiplication (Hadamard product)\n",
    "        F = np.multiply(A, B)\n",
    "        print(\"\\nElement-wise Multiplication (A .* B):\\n\", F)\n",
    "        # Inversion of matrix A\n",
    "        if np.allclose(np.linalg.det(A), 0, atol=1e-8):  # Check if A is invertible\n",
    "            print(\"\\nMatrix A is not invertible.\")\n",
    "        else:\n",
    "            A_inv = inv(A)\n",
    "            print(\"\\nInversion of A:\\n\", A_inv)\n",
    "        # Eigenvalues and eigenvectors of matrix A\n",
    "        eigenvalues, eigenvectors = eig(A)\n",
    "        print(\"\\nEigenvalues of A:\\n\", eigenvalues)\n",
    "        print(\"\\nEigenvectors of A:\\n\", eigenvectors)\n",
    "        # Singular Value Decomposition of matrix A\n",
    "        U, s, Vt = svd(A, full_matrices=False)\n",
    "        print(\"\\nSingular values of A:\\n\", s)\n",
    "        print(\"\\nLeft singular vectors of A:\\n\", U)\n",
    "        print(\"\\nRight singular vectors of A:\\n\", Vt)\n",
    "        # Selection\n",
    "        # Select the first row of matrix A\n",
    "        first_row = A[0, :]\n",
    "        print(\"\\nFirst row of A:\\n\", first_row)\n",
    "        # Select the second column of matrix A\n",
    "        second_column = A[:, 1]\n",
    "        print(\"\\nSecond column of A:\\n\", second_column)\n",
    "        # Select a submatrix from A consisting of rows 1 and 2 and columns 1 and 2\n",
    "        submatrix = A[1:3, 1:3]\n",
    "        print(\"\\nSubmatrix of A (rows 1 and 2, columns 1 and 2):\\n\", submatrix)\n",
    "\n",
    "        return 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "for i in range(5):\n",
    "    print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n"
     ]
    }
   ],
   "source": [
    "print(list(range(10)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 1.04719755, 2.0943951 , 3.14159265, 4.1887902 ,\n",
       "       5.23598776])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(0, 360, num=6, endpoint=False)/180*constv.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1., 2., 3., 4., 5.])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(0, 5, num=6, endpoint=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6, 91, 181)\n"
     ]
    }
   ],
   "source": [
    " print((x*(np.sin(phi_grid) * np.cos(theta_grid))).shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(91, 181)\n"
     ]
    }
   ],
   "source": [
    " print((np.sin(phi_grid) * np.cos(theta_grid)).shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6, 1, 1)\n"
     ]
    }
   ],
   "source": [
    "print(x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        ],\n",
       "       [0.00194805],\n",
       "       [0.0038961 ],\n",
       "       [0.00584416],\n",
       "       [0.00779221],\n",
       "       [0.00974026],\n",
       "       [0.01168831],\n",
       "       [0.01363636],\n",
       "       [0.01558442],\n",
       "       [0.01753247],\n",
       "       [0.01948052],\n",
       "       [0.02142857],\n",
       "       [0.02337662],\n",
       "       [0.02532468],\n",
       "       [0.02727273],\n",
       "       [0.02922078]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array_geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4.70097342 +4.65291933j, -1.9010936  -2.48624112j,\n",
       "        -1.6731716  +3.11916711j, ..., -4.0016807  +0.98173058j,\n",
       "         3.64336533 -4.73823158j,  8.93252885 -6.41762679j],\n",
       "       [ 7.57741845 +0.68490826j, -2.80233393 -1.4728114j ,\n",
       "        -1.78588623 +4.75358954j, ..., -3.97901989 +1.71831306j,\n",
       "         2.25326739 -3.77172293j,  4.4053539  -8.47597453j],\n",
       "       [ 7.97745863 -4.19966089j, -3.04691942 -0.13885044j,\n",
       "        -1.03418314 +6.20926249j, ..., -3.57745563 +2.33622359j,\n",
       "         1.5664352  -2.22421758j, -0.53272885 -7.88652811j],\n",
       "       [ 5.78496345 -8.5828119j , -2.56384762 +1.1283968j ,\n",
       "         0.36372045 +7.06360833j, ..., -2.91356087 +2.65608453j,\n",
       "         1.78225414 -0.54495184j, -4.44820836 -4.82040212j],\n",
       "       [ 1.63640789-11.19212879j, -1.49335285 +1.96105211j,\n",
       "         2.00201719 +7.06861296j, ..., -2.18006233 +2.58504118j,\n",
       "         2.83807261 +0.77858834j, -6.20443232 -0.16768401j],\n",
       "       [-3.26389445-11.27013446j, -0.14619664 +2.11739809j,\n",
       "         3.40511431 +6.22282354j, ..., -1.5898925  +2.14371723j,\n",
       "         4.42738958 +1.3621831j , -5.29157401 +4.72095569j]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array_signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[14724.20666228-2.62831879e-14j, 13633.03484348+3.84514707e+03j,\n",
       "        10705.1027692 +6.56599726e+03j,  6790.38060173+7.37269566e+03j,\n",
       "         3025.30079633+6.03106009e+03j,   502.85506596+2.93056343e+03j],\n",
       "       [13633.03484348-3.84514707e+03j, 14864.35425509+1.46355242e-14j,\n",
       "        13921.80718779+3.92592773e+03j, 11079.01231284+6.79295072e+03j,\n",
       "         7161.2247113 +7.76878053e+03j,  3305.76672502+6.57013662e+03j],\n",
       "       [10705.1027692 -6.56599726e+03j, 13921.80718779-3.92592773e+03j,\n",
       "        15301.75135088+2.46858816e-14j, 14444.34148292+4.07210053e+03j,\n",
       "        11598.48117391+7.10825493e+03j,  7590.31539148+8.22707641e+03j],\n",
       "       [ 6790.38060173-7.37269566e+03j, 11079.01231284-6.79295072e+03j,\n",
       "        14444.34148292-4.07210053e+03j, 15909.42284536+4.03632117e-14j,\n",
       "        15048.94754332+4.24123194e+03j, 12112.70905343+7.42037801e+03j],\n",
       "       [ 3025.30079633-6.03106009e+03j,  7161.2247113 -7.76878053e+03j,\n",
       "        11598.48117391-7.10825493e+03j, 15048.94754332-4.24123194e+03j,\n",
       "        16510.9634717 +0.00000000e+00j, 15560.10998878+4.38422361e+03j],\n",
       "       [  502.85506596-2.93056343e+03j,  3305.76672502-6.57013662e+03j,\n",
       "         7590.31539148-8.22707641e+03j, 12112.70905343-7.42037801e+03j,\n",
       "        15560.10998878-4.38422361e+03j, 16931.74773627+0.00000000e+00j]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(array_signal,(array_signal.conj().T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0, 10])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DOA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "helloworld\n",
      "[[[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849522-5.48387854e-02j\n",
      "    0.99398542-1.09512530e-01j  0.98648415-1.63856691e-01j\n",
      "    0.97601401-2.17707716e-01j  0.9626065 -2.70903537e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99398725-1.09495918e-01j\n",
      "    0.97602129-2.17675092e-01j  0.94631818-3.23236612e-01j\n",
      "    0.90523511-4.24911047e-01j  0.85326613-5.21475711e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999997+2.59781096e-04j\n",
      "    0.99999987-5.19562173e-04j -0.9999997 +7.79343216e-04j\n",
      "    0.99999946-1.03912421e-03j -0.99999916+1.29890513e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999856+1.69606750e-03j\n",
      "    0.99999425-3.39213013e-03j -0.99998706+5.08818299e-03j\n",
      "    0.99997699-6.78422122e-03j -0.99996404+8.48023993e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999763+2.17487762e-03j\n",
      "    0.99999054-4.34974496e-03j -0.99997871+6.52459172e-03j\n",
      "    0.99996216-8.69940762e-03j -0.99994087+1.08741824e-02j]]\n",
      "\n",
      " [[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849568-5.48304416e-02j\n",
      "    0.99398725-1.09495918e-01j  0.98648826-1.63831961e-01j\n",
      "    0.97602129-2.17675092e-01j  0.96261782-2.70863317e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99398907-1.09479308e-01j\n",
      "    0.97602856-2.17642473e-01j  0.94633438-3.23189172e-01j\n",
      "    0.90526351-4.24850539e-01j  0.85330969-5.21404417e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999998-2.18738244e-04j\n",
      "    0.9999999 +4.37476477e-04j -0.99999978-6.56214689e-04j\n",
      "    0.99999962+8.74952869e-04j -0.9999994 -1.09369101e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999926+1.21732992e-03j\n",
      "    0.99999704-2.43465803e-03j -0.99999333+3.65198254e-03j\n",
      "    0.99998814-4.86930163e-03j -0.99998148+6.08661351e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999856+1.69606750e-03j\n",
      "    0.99999425-3.39213013e-03j -0.99998706+5.08818299e-03j\n",
      "    0.99997699-6.78422122e-03j -0.99996404+8.48023993e-03j]]\n",
      "\n",
      " [[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849705-5.48054126e-02j\n",
      "    0.99399273-1.09446086e-01j  0.98650058-1.63757776e-01j\n",
      "    0.97604311-2.17577229e-01j  0.96265176-2.70742667e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99399456-1.09429484e-01j\n",
      "    0.97605038-2.17544623e-01j  0.94638297-3.23046861e-01j\n",
      "    0.90534867-4.24669023e-01j  0.85344035-5.21190537e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999863-1.65414976e-03j\n",
      "    0.99999453+3.30829499e-03j -0.99998769-4.96243118e-03j\n",
      "    0.99997811+6.61655378e-03j -0.9999658 -8.27065828e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999998-2.18738244e-04j\n",
      "    0.9999999 +4.37476477e-04j -0.99999978-6.56214689e-04j\n",
      "    0.99999962+8.74952869e-04j -0.9999994 -1.09369101e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999997+2.59781096e-04j\n",
      "    0.99999987-5.19562173e-04j -0.9999997 +7.79343216e-04j\n",
      "    0.99999946-1.03912421e-03j -0.99999916+1.29890513e-03j]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849705+5.48054126e-02j\n",
      "    0.99399273+1.09446086e-01j  0.98650058+1.63757776e-01j\n",
      "    0.97604311+2.17577229e-01j  0.96265176+2.70742667e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99399456+1.09429484e-01j\n",
      "    0.97605038+2.17544623e-01j  0.94638297+3.23046861e-01j\n",
      "    0.90534867+4.24669023e-01j  0.85344035+5.21190537e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999863+1.65414976e-03j\n",
      "    0.99999453-3.30829499e-03j -0.99998769+4.96243118e-03j\n",
      "    0.99997811-6.61655378e-03j -0.9999658 +8.27065828e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999998+2.18738244e-04j\n",
      "    0.9999999 -4.37476477e-04j -0.99999978+6.56214689e-04j\n",
      "    0.99999962-8.74952869e-04j -0.9999994 +1.09369101e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999997-2.59781096e-04j\n",
      "    0.99999987+5.19562173e-04j -0.9999997 -7.79343216e-04j\n",
      "    0.99999946+1.03912421e-03j -0.99999916-1.29890513e-03j]]\n",
      "\n",
      " [[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849568+5.48304416e-02j\n",
      "    0.99398725+1.09495918e-01j  0.98648826+1.63831961e-01j\n",
      "    0.97602129+2.17675092e-01j  0.96261782+2.70863317e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99398907+1.09479308e-01j\n",
      "    0.97602856+2.17642473e-01j  0.94633438+3.23189172e-01j\n",
      "    0.90526351+4.24850539e-01j  0.85330969+5.21404417e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999998+2.18738244e-04j\n",
      "    0.9999999 -4.37476477e-04j -0.99999978+6.56214689e-04j\n",
      "    0.99999962-8.74952869e-04j -0.9999994 +1.09369101e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999926-1.21732992e-03j\n",
      "    0.99999704+2.43465803e-03j -0.99999333-3.65198254e-03j\n",
      "    0.99998814+4.86930163e-03j -0.99998148-6.08661351e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999856-1.69606750e-03j\n",
      "    0.99999425+3.39213013e-03j -0.99998706-5.08818299e-03j\n",
      "    0.99997699+6.78422122e-03j -0.99996404-8.48023993e-03j]]\n",
      "\n",
      " [[ 1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j\n",
      "    1.        +0.00000000e+00j  1.        +0.00000000e+00j]\n",
      "  [ 1.        +0.00000000e+00j  0.99849522+5.48387854e-02j\n",
      "    0.99398542+1.09512530e-01j  0.98648415+1.63856691e-01j\n",
      "    0.97601401+2.17707716e-01j  0.9626065 +2.70903537e-01j]\n",
      "  [ 1.        +0.00000000e+00j  0.99398725+1.09495918e-01j\n",
      "    0.97602129+2.17675092e-01j  0.94631818+3.23236612e-01j\n",
      "    0.90523511+4.24911047e-01j  0.85326613+5.21475711e-01j]\n",
      "  ...\n",
      "  [ 1.        +0.00000000e+00j -0.99999997-2.59781096e-04j\n",
      "    0.99999987+5.19562173e-04j -0.9999997 -7.79343216e-04j\n",
      "    0.99999946+1.03912421e-03j -0.99999916-1.29890513e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999856-1.69606750e-03j\n",
      "    0.99999425+3.39213013e-03j -0.99998706-5.08818299e-03j\n",
      "    0.99997699+6.78422122e-03j -0.99996404-8.48023993e-03j]\n",
      "  [ 1.        +0.00000000e+00j -0.99999763-2.17487762e-03j\n",
      "    0.99999054+4.34974496e-03j -0.99997871-6.52459172e-03j\n",
      "    0.99996216+8.69940762e-03j -0.99994087-1.08741824e-02j]]]\n",
      "IP: 127.0.0.1;                 PORT: 12345;                 POWER: 4;                 ELEMENTCNT: 6;                 ARRAYTYPE: ULA;                 SPACING: 6.25;                 FREQUENCE: 2.4;                 \n",
      "array_signal_shape: (6, 400)\n",
      "doaest\n",
      "\n",
      "Inversion of A:\n",
      " [[ 0.96191311-1.06415872e-18j -0.0429495 -1.19164685e-02j\n",
      "  -0.04092318-2.46263570e-02j -0.03259616-3.44400218e-02j\n",
      "  -0.02038576-3.85085846e-02j -0.00783661-3.56509534e-02j]\n",
      " [-0.0429495 +1.19164685e-02j  0.94769732+1.07368445e-18j\n",
      "  -0.05411512-1.50399198e-02j -0.04786069-2.88372459e-02j\n",
      "  -0.03535502-3.73866547e-02j -0.02022848-3.82062777e-02j]\n",
      " [-0.04092318+2.46263570e-02j -0.05411512+1.50399198e-02j\n",
      "   0.93958187+5.83348482e-18j -0.05800246-1.61273574e-02j\n",
      "  -0.04756938-2.86604275e-02j -0.03214757-3.39608965e-02j]\n",
      " [-0.03259616+3.44400218e-02j -0.04786069+2.88372459e-02j\n",
      "  -0.05800246+1.61273574e-02j  0.93992264+3.87230323e-18j\n",
      "  -0.05348304-1.48631014e-02j -0.04013382-2.41472317e-02j]\n",
      " [-0.02038576+3.85085846e-02j -0.03535502+3.73866547e-02j\n",
      "  -0.04756938+2.86604275e-02j -0.05348304+1.48631014e-02j\n",
      "   0.94862072+1.22934006e-18j -0.04186882-1.16141616e-02j]\n",
      " [-0.00783661+3.56509534e-02j -0.02022848+3.82062777e-02j\n",
      "  -0.03214757+3.39608965e-02j -0.04013382+2.41472317e-02j\n",
      "  -0.04186882+1.16141616e-02j  0.96315107+4.52596938e-18j]]\n",
      "U\n",
      " [[ 0.96191311-1.06415872e-18j -0.0429495 -1.19164685e-02j\n",
      "  -0.04092318-2.46263570e-02j -0.03259616-3.44400218e-02j\n",
      "  -0.02038576-3.85085846e-02j -0.00783661-3.56509534e-02j]\n",
      " [-0.0429495 +1.19164685e-02j  0.94769732+1.07368445e-18j\n",
      "  -0.05411512-1.50399198e-02j -0.04786069-2.88372459e-02j\n",
      "  -0.03535502-3.73866547e-02j -0.02022848-3.82062777e-02j]\n",
      " [-0.04092318+2.46263570e-02j -0.05411512+1.50399198e-02j\n",
      "   0.93958187+5.83348482e-18j -0.05800246-1.61273574e-02j\n",
      "  -0.04756938-2.86604275e-02j -0.03214757-3.39608965e-02j]\n",
      " [-0.03259616+3.44400218e-02j -0.04786069+2.88372459e-02j\n",
      "  -0.05800246+1.61273574e-02j  0.93992264+3.87230323e-18j\n",
      "  -0.05348304-1.48631014e-02j -0.04013382-2.41472317e-02j]\n",
      " [-0.02038576+3.85085846e-02j -0.03535502+3.73866547e-02j\n",
      "  -0.04756938+2.86604275e-02j -0.05348304+1.48631014e-02j\n",
      "   0.94862072+1.22934006e-18j -0.04186882-1.16141616e-02j]\n",
      " [-0.00783661+3.56509534e-02j -0.02022848+3.82062777e-02j\n",
      "  -0.03214757+3.39608965e-02j -0.04013382+2.41472317e-02j\n",
      "  -0.04186882+1.16141616e-02j  0.96315107+4.52596938e-18j]]\n",
      "pu.shape (181, 91)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gl/venv0825/lib/python3.6/site-packages/numpy/core/_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/ipykernel_launcher.py:143: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/ipykernel_launcher.py:147: RuntimeWarning: divide by zero encountered in log10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoJUlEQVR4nO3deXxU9f398dc7OyEJEBL2kIWEVRAwArIIVquIKJWqBfddq6ittVa/rXXpotVarYq1WHfrWq2g4oICFSsqYV8DIQImQBL2sCQhyef3R4b+xhggwCR3MjnPx2MemfncO5PDTThzc++de805h4iIhK4wrwOIiEjDUtGLiIQ4Fb2ISIhT0YuIhDgVvYhIiFPRi4iEuMMWvZk9a2bFZrbsINPNzB4zszwzW2JmA/2mVZnZIt9tWiCDi4hI/dRnjf55YPQhpp8JZPlu1wJ/85u2zznX33c756hTiojIUTts0TvnPgO2HWKWccCLrsaXQGsz6xiogCIicmwiAvAanYFv/R4X+MY2ATFmlgNUAg8459453IslJSW5tLS0AMQSEWk+5s+fv8U5l1zXtEAU/aGkOucKzSwDmGlmS51za2vPZGbXUrPZh65du5KTk9PAsUREQouZrT/YtEAcdVMIpPg97uIbwzl34Gs+MBsYUNcLOOemOOeynXPZycl1viGJiMhRCkTRTwMu9R19MwTY6ZzbZGZtzCwawMySgGHAigB8PxEROQKH3XRjZq8Co4AkMysA7gYiAZxzTwHTgTFAHrAXuML31F7A382smpo3lAeccyp6EZFGdtiid85NPMx0B9xYx/gXQN+jjyYiIoGgT8aKiIQ4Fb2ISIhT0YuIhLiGPo6+0eytqORvs9diZoSbER5Gzf2wmsdm1NwPM8yM6PAwoiPDiI4IJzoyjJhaX2OjwomPiaRlVDhm5vU/T0TkqIVQ0VcxeVYe1QG+BG6YQXxMJPExEST4vraJjaJdQjTt4qNJjo+mXXxMzdeEaJJaRhMWpjcGEQkeIVP0SXHR5N9/Fs45qh1UO0dVtaPa97iq2lHte1zlHBWV1ZRXVlO+v5qyyqrvfC2vrGJfRRWlZZXsKtvPrn37/e5Xkleym7n5W9m5b//3ckRHhJGSGEtqYiypbVuS2jaWzHZx9OqYQGLLKA+WjIg0dyFT9AfUbLqBcIzI8Ib9XmX7qygpLadkdznFu8opLi2jYPs+1m/dw/qte/li7Vb27a/63/ztE6Lp1TGBXh0TOL5LK7LTEkmKi27YkCLS7IVc0TemmMhwUhJjSUmMrXO6c46S0nJWF+1m5aZdrNy0ixWbdvH5mi1U+rYxZSS3ZHB6IiemJXJSt7Z0bNWiMf8JItIMqOgbkJnRLiGGdgkxDM9K+t94eWUVywp38vU325m3bhvvLdnEq1/XnAC0V8cEftAzmR/0bM+AlNba3i8ix8xqPtgaPLKzs11zO3tldbVj1eZS5qwpYeaqYnLWb6eq2tGxVQxj+3VkbL9O9OvSSkf/iMhBmdl851x2ndNU9MFn5979zMot5r0lG/nP6hL2Vzky28Ux4cQUfjywC220U1dEalHRN2E79+7ng2WbeD3nWxZu2EFUeBhn9evI1SPS6dOpldfxRCRIqOhDxKrNu3j1qw38a34BeyqqGJ6ZxLUnZzAiK0mbdUSaORV9iNm5bz+vfLWB5/77DcWl5Qzs2ppfntGTk7q19TqaiHhERR+iyiureGt+IY99uobNu8oYnpnEr0b3pG8XbdIRaW4OVfQ6qVkTFh0RzoWDuzL7l6O4a2xvVm7axTmTP+fOt5ewbU+F1/FEJEio6ENATGQ4Vw1PZ/YvR3HVsHTezCnglD/P5qW566gO9Ml/RKTJUdGHkPiYSH4ztjcf3DKC4zoncNfU5Ux4+kvWb93jdTQR8ZCKPgRltY/n5asG89B5/Vi5aRejH53Dc//9Rmv3Is2Uij5EmRnnZ6cw4+cjGZKRyL3vruCy575my+5yr6OJSCNT0Ye4Dq1iePbyE/njuX356pttjPnrHL7M3+p1LBFpRCr6ZsDMuHBwV965YRhx0RFc+PSXTJ6VR7AdWisiDUNF34z07pTAtJuGM7ZfJx76KJebX1tEmd/58kUkNOk0xc1MXHQEf53Qn54d43noo1w2bN3DlEuzaZ8Q43U0EWkgWqNvhsyMG0Zl8veLT2BN8W7OeeJzVm7a5XUsEWkghy16M3vWzIrNbNlBppuZPWZmeWa2xMwG+k27zMzW+G6XBTK4HLvT+3TgrZ8OxTB+8ve55Kzb5nUkEWkA9Vmjfx4YfYjpZwJZvtu1wN8AzCwRuBsYDAwC7jazNscSVgKvV8cE/vXTk0iKi+biZ75i1qpiryOJSIAdtuidc58Bh1rVGwe86Gp8CbQ2s47AGcAM59w259x2YAaHfsMQj3RpE8sb159EZrs4rnkxh3cXb/Q6kogEUCC20XcGvvV7XOAbO9j495jZtWaWY2Y5JSUlAYgkRyopLppXrxnCwNQ23PLaQt5borIXCRVBsTPWOTfFOZftnMtOTk72Ok6zFR8TyXOXn0h2aiK3vLaI6Us3eR1JRAIgEEVfCKT4Pe7iGzvYuASxltERPHvFiQxIac3Nry7kw2WbvY4kIscoEEU/DbjUd/TNEGCnc24T8BFwupm18e2EPd03JkEuLjqC5644kX5dWnHTqwv4fM0WryOJyDGoz+GVrwJzgR5mVmBmV5nZ9WZ2vW+W6UA+kAc8DdwA4JzbBvwOmOe73ecbkyYgPiaS564YRLfkOK57KYdlhTu9jiQiR0mXEpRDKtpVxvgnv6i5bOFPh5LatqXXkUSkDrqUoBy19gkxvHDlIKqqHZc88zUlpTrNsUhTo6KXw8psF8czl59IcWkZ1788n/JKnQhNpClR0Uu9DOzahofP78/89du58+2lOsWxSBOiopd6O6tfR352WhZvLyhkymf5XscRkXrSaYrliNz8gyzWFO3mgQ9X0S05jtN6t/c6kogchtbo5YiEhRl/Pv94juvUip+/voh1W/Z4HUlEDkNFL0esRVQ4f7t4IOHhxk//uUBXqRIJcip6OSpd2sTyyE/6s2rzLu56p85LFYhIkFDRy1E7pUc7bjolkzfnF/D6vA1exxGRg1DRyzG55bTujMhK4q6py1m1WZcjFAlGKno5JuFhxiM/6U9CTCQ3v7pQ2+tFgpCKXo5ZUlw0D19wPKuLdnP/9JVexxGRWlT0EhAjuydz1fB0Xpi7npmriryOIyJ+VPQSMLeP7kGvjgn88s0lFJeWeR1HRHxU9BIw0RHhPDahP7vLK7nzLZ0PRyRYqOgloLLax3P76J58uqqYtxfoypEiwUBFLwF3xdA0BqUlcs+7y9m8U5twRLymopeACwszHjyvH5VVjjveXqJNOCIeU9FLg0hLaskdZ/Zkdm4Jb+YUeB1HpFlT0UuDuWRIKkMyEvndeyso2qVNOCJeUdFLgwkLMx4Y34+KqmrufXe513FEmi0VvTSotKSW3HxqFtOXbuaTFfoglYgXVPTS4K4ZkUGP9vH8duoydpdXeh1HpNlR0UuDi4oI44/j+7JpVxkPf5zrdRyRZkdFL43ihNQ2XDIklee/WMeSgh1exxFpVupV9GY22sxyzSzPzO6oY3qqmX1qZkvMbLaZdfGbVmVmi3y3aYEML03LL8/oQVJcNHdNXU51tY6tF2kshy16MwsHJgNnAr2BiWbWu9ZsfwZedM71A+4D7vebts851993OydAuaUJio+J5P/G9GTxtzt4I+dbr+OINBv1WaMfBOQ55/KdcxXAa8C4WvP0Bmb67s+qY7oIAD/q35lBaYn86cNV7Nhb4XUckWahPkXfGfBf/SrwjflbDIz33T8XiDeztr7HMWaWY2ZfmtmP6voGZnatb56ckpKS+qeXJsfMuHdcH3aVVfJn7ZgVaRSB2hl7GzDSzBYCI4FC4MA15VKdc9nAhcCjZtat9pOdc1Occ9nOuezk5OQARZJg1atjApcMSeWfX21gacFOr+OIhLz6FH0hkOL3uItv7H+ccxudc+OdcwOAX/vGdvi+Fvq+5gOzgQHHnFqavFtP707bltHc8+5ynfRMpIHVp+jnAVlmlm5mUcAE4DtHz5hZkpkdeK07gWd9423MLPrAPMAwYEWgwkvTlRATyW2nd2f++u1MX7rZ6zgiIe2wRe+cqwQmAR8BK4E3nHPLzew+MztwFM0oINfMVgPtgT/4xnsBOWa2mJqdtA8451T0AsD52Sn07BDP/R+spGx/1eGfICJHxYLtz+bs7GyXk5PjdQxpJJ+v2cLFz3zFHWf25PqR39t9IyL1ZGbzfftDv0efjBVPDc9K4tSe7Zg8M48tu8u9jiMSklT04rn/O6sX+/ZX8ciM1V5HEQlJKnrxXLfkOC4eksqrX29gdVGp13FEQo6KXoLCLadmERcdwe/fX+l1FJGQo6KXoNCmZRQ3n5rFZ6tLmJ1b7HUckZCiopegcelJaaS1jeX+6auo0tktRQJGRS9BIyoijNvO6EFuUSlTFxUe/gkiUi8qegkqY47ryHGdE/jLjNWUV+pDVCKBoKKXoBIWZtx+Rk8Ktu/j1a82eB1HJCSo6CXojMhK4qSMtjw+M08XExcJABW9BB0z41dn9mTrngqemfON13FEmjwVvQSl/imtGd2nA0/PyWerTo0gckxU9BK0bjujO3srKnly9lqvo4g0aSp6CVqZ7eI574QuvDR3PYU79nkdR6TJUtFLULvltO5g6IRnIsdARS9BrXPrFlwyJJW3FxTwzZY9XscRaZJU9BL0rh/ZjaiIMB77dI3XUUSaJBW9BL3k+GguOymNqYsKySvWaYxFjpSKXpqEa0/OICYynL9+mud1FJEmR0UvTULbuGguH5rGe0s2krtZa/UiR0JFL03GNSMyaBkVwV8/1RE4IkdCRS9NRpuWUVw5LI3pSzezYuMur+OINBkqemlSrhqeQXxMBI9+orV6kfpS0UuT0io2kquHZ/DxiiKWFuz0Oo5Ik1Cvojez0WaWa2Z5ZnZHHdNTzexTM1tiZrPNrIvftMvMbI3vdlkgw0vzdMXwNFq1iOQRrdWL1Mthi97MwoHJwJlAb2CimfWuNdufgRedc/2A+4D7fc9NBO4GBgODgLvNrE3g4ktzlBATyTUj0pm5qlhr9SL1UJ81+kFAnnMu3zlXAbwGjKs1T29gpu/+LL/pZwAznHPbnHPbgRnA6GOPLc3dpUPTSIiJ4PGZ+rSsyOHUp+g7A9/6PS7wjflbDIz33T8XiDeztvV8rsgRS4iJ5Iph6Xy8ooiVm3QEjsihBGpn7G3ASDNbCIwECoF6X9nZzK41sxwzyykpKQlQJAl1Vw5LJy46gidm6tOyIodSn6IvBFL8Hnfxjf2Pc26jc268c24A8Gvf2I76PNc37xTnXLZzLjs5OfnI/gXSbLWKjeSyoalMX7aJNUX6tKzIwdSn6OcBWWaWbmZRwARgmv8MZpZkZgde607gWd/9j4DTzayNbyfs6b4xkYC4angGLSLDeWKW1upFDuawRe+cqwQmUVPQK4E3nHPLzew+MzvHN9soINfMVgPtgT/4nrsN+B01bxbzgPt8YyIBkdgyikuGpPLu4o3kl+z2Oo5IUDLnnNcZviM7O9vl5OR4HUOakJLSckY8OJOz+nbi4QuO9zqOiCfMbL5zLruuafpkrDR5yfHRXDgolXcWFbJh616v44gEHRW9hITrRmYQHmY8OVvb6kVqU9FLSGifEMOEE1N4a0EBBdu1Vi/iT0UvIeP6kd0AeOo/az1OIhJcVPQSMjq1bsF5J6TwxrwCNu8s8zqOSNBQ0UtIuWFUN6qd01q9iB8VvYSUlMRYzh3QmVe/3kBJabnXcUSCgopeQs5PR3Vjf1U1//g83+soIkFBRS8hJyM5jrH9OvHy3PVs31PhdRwRz6noJSRN+kEmeyqqeO6/33gdRcRzKnoJSd3bxzO6Twee+2Idu8r2ex1HxFMqeglZk36QSWlZJS/NXe91FBFPqeglZB3XuRWn9EjmH3Py2VtR6XUcEc+o6CWkTfpBFtv37uefX27wOoqIZ1T0EtJOSG3DsMy2TJmTT9n+el/dUiSkqOgl5E06JYuS0nLeyPn28DOLhCAVvYS8IRmJZKe24anZa6morPY6jkijU9FLyDMzbjo1i407y3h7QYHXcUQanYpemoWTs5Lo16UVT85eS2WV1uqleVHRS7NgZkw6JZMN2/by7pKNXscRaVQqemk2TuvVnp4d4nliZh7V1c7rOCKNRkUvzUZYmDHpB5msLdnDB8s2ex1HpNGo6KVZOfO4jmQkt+TxmWtwTmv10jyo6KVZCQ8zbhyVyarNpXy6stjrOCKNQkUvzc45/TuRkthCa/XSbNSr6M1stJnlmlmemd1Rx/SuZjbLzBaa2RIzG+MbTzOzfWa2yHd7KtD/AJEjFRkexg2jMllcsJM5a7Z4HUekwR226M0sHJgMnAn0BiaaWe9as/0GeMM5NwCYADzpN22tc66/73Z9gHKLHJPxAzvTsVUMT8zM8zqKSIOrzxr9ICDPOZfvnKsAXgPG1ZrHAQm++60AHagsQS06IpzrTs7g63Xb+Cp/q9dxRBpUfYq+M+B/NqgC35i/e4CLzawAmA7c5Dct3bdJ5z9mNqKub2Bm15pZjpnllJSU1D+9yDGYMKgrSXHRPK61eglxgdoZOxF43jnXBRgDvGRmYcAmoKtvk86twCtmllD7yc65Kc65bOdcdnJycoAiiRxaTGQ4156czud5W1iwYbvXcUQaTH2KvhBI8XvcxTfm7yrgDQDn3FwgBkhyzpU757b6xucDa4HuxxpaJFAuGpxK69hIJmutXkJYfYp+HpBlZulmFkXNztZptebZAJwKYGa9qCn6EjNL9u3MxcwygCwgP1DhRY5Vy+gIrhqWzqerillWuNPrOCIN4rBF75yrBCYBHwErqTm6ZrmZ3Wdm5/hm+wVwjZktBl4FLnc1ByifDCwxs0XAv4DrnXPbGuDfIXLULh2aRnx0BJNnaa1eQlNEfWZyzk2nZier/9hv/e6vAIbV8by3gLeOMaNIg2rVIpLLh6Xx+Mw8VheV0r19vNeRRAJKn4wVAa4Ylk5sVDhPaq1eQpCKXgRIbBnFxUNSmbZ4I+u27PE6jkhAqehFfK4ekU5keBhPztZavYQWFb2IT7v4GCYO6spbCwpZv1Vr9RI6VPQifm4Y1Y2IMNOnZSWkqOhF/LRLiOGSIam8vaCA/JLdXscRCQgVvUgt143sRnREOI99usbrKCIBoaIXqSU5PppLh6YydfFG8opLvY4jcsxU9CJ1uO7kbsRGhvPIJ1qrl6ZPRS9Sh8SWUVwxLJ33l2xi1eZdXscROSYqepGDuHpEOvHRETw6Q2v10rSp6EUOonVsFFcOT+fD5Zt1Zktp0lT0Iodw1Yh0EmIiePST1V5HETlqKnqRQ0iIieS6kd34ZGUx89frDNvSNNXrNMUizdkVw9J4/ot1/OmDXF6/bghm5nUkaYK27C5n/dY9bNi2lw1b97Fldzk79+1n57797C6vpLKqmsx28Tx8wfEB/94qepHDiI2K4OZTs7jrnWXMXFXMqb3aex1JgtzeikrmrdvO/HXbWLZxF0sLd1JSWv6dedrERtKqRc0tLiaCyJgI2sRGNkgeFb1IPUw4MYVn5uTz4Ie5jOrRjvAwrdXLd+UVl/Lhss18tmYLCzdsZ3+VI8wgs10cIzKT6NO5FRnJLemaGEvn1i2IiQxvtGwqepF6iAwP47YzejDplYW8s7CQH5/QxetIEgTyS3bzzsJCpi/bTF5xzbmRjuucwJXD0jmpW1tOTEukZbT3Net9ApEmYsxxHenbOZ+/zFjNWf06NuoamQSPvRWVvL9kE2/mFPD1um2EGQxKT+SSIX04o08HOrSK8Tri96joReopLMz41eieXPzMV7z85XquHpHhdSRpRJt27uOFL9bz6tcb2LlvPxlJLfnV6J6MH9iZ9gnBV+7+VPQiR2B4VhLDM5OYPCuPC05MISGmYXaeSfDIKy5l8qy1vLt4I9XOcUafDlw+NI1B6YlN5ggsFb3IEfrV6J6c/cTnPP1ZPr84vYfXcaSB5BWX8tineby7ZCMtIsO59KQ0rhiWRkpirNfRjpiKXuQI9e3SirH9OvKPOd9w0eDUoNwmK0evaFcZD32Uy1sLCmgRGc71I7txzYgMEltGeR3tqKnoRY7C7Wf05OPlRTz0UW6DfMBFGt++iiqmfJbPU/9ZS1W14+rh6fx0VGaTLvgDVPQiR6Fr21iuGJ7G3/+Tz2VDU+nXpbXXkeQoVVc7pi4u5MEPc9m0s4wxfTtwx+hedG3b9DbRHEy9znVjZqPNLNfM8szsjjqmdzWzWWa20MyWmNkYv2l3+p6Xa2ZnBDK8iJcmnZJJUlwUv3tvBc45r+PIUVhTVMoFf5/Lz19fTHJ8NG9cdxJPXnRCSJU81KPozSwcmAycCfQGJppZ71qz/QZ4wzk3AJgAPOl7bm/f4z7AaOBJ3+uJNHnxMZHc+sMezFu3nelLN3sdR45AeWUVj8xYzZjH5pBXspsHz+vHOzcMY1B6otfRGkR91ugHAXnOuXznXAXwGjCu1jwOSPDdbwVs9N0fB7zmnCt3zn0D5PleTyQk/OTEFHp2iOf+D1ZStr/K6zhSD/PWbeOsxz7nr5+u4ay+Hfnk1pFckJ1CWAif1qI+Rd8Z+NbvcYFvzN89wMVmVgBMB246gudiZteaWY6Z5ZSUlNQzuoj3wsOMu8b2pmD7Pp7+LN/rOHII+yqquHvqMs5/ai5l+6t4/ooTeXTCAJLior2O1uACdT76icDzzrkuwBjgJTOr92s756Y457Kdc9nJyckBiiTSOIZlJnHmcR2YPDuPb7ft9TqO1GFJwQ7OenwOL8xdz5XD0vn45yczqkc7r2M1mvqUcSGQ4ve4i2/M31XAGwDOublADJBUz+eKNHl3je2NYdz33gqvo4ifyqpqHvt0DeOf/IKyiipeuXowvz27N7FRzeuAw/oU/Twgy8zSzSyKmp2r02rNswE4FcDMelFT9CW++SaYWbSZpQNZwNeBCi8SLDq1bsHNp2YxY0URM1cVeR1HgHVb9nDeU3P5y4zVjO3XkQ9+djJDM5O8juWJw76tOecqzWwS8BEQDjzrnFtuZvcBOc65acAvgKfN7OfU7Ji93NUcb7bczN4AVgCVwI3OOe2xkpB01fB03pz/LfdMW8HQbkk6u6WHpi4q5P/eXkpEeBiPTxzA2cd38jqSpyzYjv/Nzs52OTk5XscQOSr/zdvCRf/4ip+dlsXPTuvudZxmZ19FFfe+u5zX5n1LdmobHps4gE6tW3gdq1GY2XznXHZd05rXhiqRBjYsM4mzj+/Ek7PWMrZfRzLbxXsdqdnIKy7lxn8uJLeolBtGdePWH3YnIjxQx5s0bVoKIgH227G9aREVzh1vLaW6Orj+Yg5Vb80v4OzH/8uW3eW8cOUgbh/dUyXvR0tCJMCS46O5a2xvctZv5+Wv1nsdJ6RVVFbzm3eW8os3F3N8Sis+uGUEI7vrEO3aVPQiDeDHAzszIiuJP32wisId+7yOE5I27yzjJ1Pm8vKXG7ju5Axevmow7YL8Sk9eUdGLNAAz44/n9qXawa//vVQnPQuwL/O3MvbxOazeXMqTFw3kzjG9tKnmELRkRBpISmIsvzyjB7NzS3gzp8DrOCHBOcc/5uRz0T++IiEmknduHMaYvh29jhX0VPQiDejyoWkMyUjk3neXs2GrTo9wLPZWVHLza4v4/fsrOa1XO6ZOGkZWex3VVB8qepEGFBZmPHxBf8LCjJ+/sYjKqmqvIzVJ32zZw7mTv+D9JRu5fXQPnrr4BOJ1YfZ6U9GLNLDOrVvwu3HHMX/9dp76z1qv4zQ5n6wo4pwnPqe4tIwXrhzEDaMyMQvdUwo3BBW9SCMY178TY/t15NFP1rCkYIfXcZqEqmrHXz7O5eoXc0htG8u7Nw1nRJYOnTwaKnqRRmBm/OFHfUmOj+bGVxawc+9+ryMFtR17K7jy+Xk8NjOP80/owr+uH0qXNqF1eb/GpKIXaSStYiN54sKBbNpRxi/eXKxPzR7E8o07OfuJz/li7Rb+cO5xPHheP50g7hip6EUa0Qmpbfj1Wb34ZGURU+boilS1/XthAeOf/IL9lY7XrzuJiwanant8AOikZiKN7PKhaeSs386DH66if0prhmS09TqS5yoqq/nD+yt4Ye56Bqcn8sSFA0mOD/1L/DUWrdGLNDIz408/7kdaUksmvbKg2V9+sGhXGROf/pIX5q7n6uHpvHz1YJV8gKnoRTwQFx3BlEuyqais5srn57GrrHnunJ23bhtjH/+cFRt38fjEAfxmbG8idSqDgNMSFfFIZrs4nrr4BL7Zsocb/7mA/c3ow1QHTmUwccqXtIwK550bhzX7q0A1JBW9iIeGZibxx/F9mbNmC7+duqxZnPxsx94Krnkxh9+/v5JTerZj6qTh9OigUxk0JO2MFfHYBdkprNuyhydnryU5LppbT+/hdaQGM3/9Nm56ZSElu8u5++zeXD40TUfVNAIVvUgQuO30HmzdXcFjM/OIjgznxlMyvY4UUNXVjilz8nnoo1w6t27BWz8dSr8urb2O1Wyo6EWCQFiY8cfxfamoquahj3KJCg/jmpMzvI4VEBt37OOX/1rMf/O2clbfjtz/474k6IRkjUpFLxIkwsOMh87rV3NM+fSVmMHVI5pu2Tvn+PfCQu6etpyqasf94/sy4cQUbarxgIpeJIhEhIfx6IT+VDvH799fSXFpOXeM7klYWNMqx627y/m/fy/lo+VFnJjWhofP70/XtjpXjVdU9CJBJjI8jCcuHMi97y5nymf5bNpZxp/P70d0RPCf78U5xwfLNvPbqcvYta+SO8/sydUjMghvYm9UoaZeRW9mo4G/AuHAP5xzD9Sa/ghwiu9hLNDOOdfaN60KWOqbtsE5d04AcouEtPAw495z+tCxVQv+9OEqineVBf1pAb7dtpffTl3GrNwS+nRK4OWrj6dnhwSvYwlghztu18zCgdXAD4ECYB4w0Tm34iDz3wQMcM5d6Xu82zkXV99A2dnZLicnp76zi4S8qYsKuf1fS4iPieTRn/RneFaS15G+o2x/Fc98/g2Pz1xDmBm3/rA7lw9N08W6G5mZzXfOZdc1rT5r9IOAPOdcvu/FXgPGAXUWPTARuPtogorI943r35keHeKZ9MpCLnn2K24clcnPTsvyvEirqx3TFm/koY9yKdyxj9N7t+eec/rQqXULT3PJ99Wn6DsD3/o9LgAG1zWjmaUC6cBMv+EYM8sBKoEHnHPvHF1UkearZ4cEpk0axt1Tl/PErDw+WVnEvef0YbAHZ750zjErt5hHZqxhaeFO+nRK4KHz+jE0M7j+0pD/L9A7YycA/3LOVfmNpTrnCs0sA5hpZkudc9+5cKaZXQtcC9C1a9cARxIJDbFRETx0/vGc2qs9v3tvBT+Z8iU/6t+JO8f0on1CTIN//6pqxwfLNjF51lpWbtpF59Yt+MsFx/Oj/p2b3FFBzU19ir4QSPF73MU3VpcJwI3+A865Qt/XfDObDQwA1taaZwowBWq20dcnuEhzNfq4DozsnszkWXlM+Syf6cs2M35AZ64ekU5mu8CfM2bTzn28mVPA6/O+pXDHPjKSW/Ln849nXP9OOtNkE1GfnbER1OyMPZWagp8HXOicW15rvp7Ah0C6872ombUB9jrnys0sCZgLjDvYjlzQzliRI7F+6x7+/lk+b80voLyymlN6JDOuf2dO6dGOVrFH/+nTjTv2MXNVMTNWFDFnTQnVDoZnJnHR4K6c3qeDDpcMQse0M9Y5V2lmk4CPqDm88lnn3HIzuw/Icc5N8806AXjNffedoxfwdzOrpuZMmQ8cquRF5Miktm3JH8/ty60/7M5Lc9fz6tcbmJVbQkSYMSSjLYPSE+nZIZ6eHRLo0qbF9zaxOOcoLa/km5I9LN+4i+UbdzJ//XZWbS4FICWxBTeMyuSC7BR94KkJO+wafWPTGr3I0auudiwu2MHHK4qYsaKIvOLd/5sWHmbERoXTMiqCFlHhlJbtZ8fe/VT6XaQ8PiaC4zq1YlSPZE7t1Y5uyXE6ZUETcag1ehW9SAjbU15JblEpuZtLKdi+lz3lVeytqGTf/mrioiNoExtJm9goUhJb0KdTK7q0aaFib6KO9Th6EWmiWkZHMLBrGwZ2beN1FPGQdpmLiIQ4Fb2ISIhT0YuIhDgVvYhIiFPRi4iEOBW9iEiIU9GLiIQ4Fb2ISIgLuk/GmlkJsP4YXiIJ2BKgOA1JOQOrqeSEppNVOQOvIbOmOueS65oQdEV/rMws52AfAw4myhlYTSUnNJ2syhl4XmXVphsRkRCnohcRCXGhWPRTvA5QT8oZWE0lJzSdrMoZeJ5kDblt9CIi8l2huEYvIiJ+QqbozWy0meWaWZ6Z3eF1ngPMLMXMZpnZCjNbbma3+MbvMbNCM1vku43xOiuAma0zs6W+TDm+sUQzm2Fma3xfPT25uZn18Ftui8xsl5n9LBiWqZk9a2bFZrbMb6zO5Wc1HvP9zi4xs4FBkPUhM1vly/NvM2vtG08zs31+y/Ypj3Me9GdtZnf6lmmumZ3hcc7X/TKuM7NFvvHGXZ7OuSZ/o+ZatmuBDCAKWAz09jqXL1tHYKDvfjw1F1rvDdwD3OZ1vjryrgOSao09CNzhu38H8Cevc9b62W8GUoNhmQInAwOBZYdbfsAY4APAgCHAV0GQ9XQgwnf/T35Z0/znC4Kcdf6sff+3FgPRQLqvF8K9yllr+sPAb71YnqGyRj8IyHPO5TvnKoDXgHEeZwLAObfJObfAd78UWAl09jbVERsHvOC7/wLwI++ifM+pwFrn3LF8yC5gnHOfAdtqDR9s+Y0DXnQ1vgRam1nHRglK3Vmdcx875yp9D78EujRWnoM5yDI9mHHAa865cufcN0AeNf3Q4A6V02quz3gB8GpjZKktVIq+M/Ct3+MCgrBMzSwNGAB85Rua5PsT+VmvN4f4ccDHZjbfzK71jbV3zm3y3d8MtPcmWp0m8N3/PMG4TA+2/IL99/ZKav7iOCDdzBaa2X/MbIRXofzU9bMO1mU6Aihyzq3xG2u05RkqRR/0zCwOeAv4mXNuF/A3oBvQH9hEzZ91wWC4c24gcCZwo5md7D/R1fzdGRSHaplZFHAO8KZvKFiX6f8E0/I7FDP7NVAJ/NM3tAno6pwbANwKvGJmCV7lown8rGuZyHdXSBp1eYZK0RcCKX6Pu/jGgoKZRVJT8v90zr0N4Jwrcs5VOeeqgadppD8vD8c5V+j7Wgz8m5pcRQc2Kfi+FnuX8DvOBBY454ogeJcpB19+Qfl7a2aXA2OBi3xvTPg2hWz13Z9Pzbbv7l5lPMTPOuiWqZlFAOOB1w+MNfbyDJWinwdkmVm6by1vAjDN40zA/7bNPQOsdM79xW/cf1vsucCy2s9tbGbW0sziD9ynZsfcMmqW5WW+2S4DpnqT8Hu+s5YUjMvU52DLbxpwqe/omyHATr9NPJ4ws9HA7cA5zrm9fuPJZhbuu58BZAH53qQ85M96GjDBzKLNLJ2anF83dr5aTgNWOecKDgw0+vJsrL2+DX2j5giG1dS8M/7a6zx+uYZT86f6EmCR7zYGeAlY6hufBnQMgqwZ1ByxsBhYfmA5Am2BT4E1wCdAYhBkbQlsBVr5jXm+TKl549kE7Kdm+/BVB1t+1BxtM9n3O7sUyA6CrHnUbOM+8Lv6lG/eH/t+JxYBC4CzPc550J818GvfMs0FzvQyp2/8eeD6WvM26vLUJ2NFREJcqGy6ERGRg1DRi4iEOBW9iEiIU9GLiIQ4Fb2ISIhT0YuIhDgVvYhIiFPRi4iEuP8HdCOOuIB3OFoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['SimSun'] not found. Falling back to DejaVu Sans.\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 26469 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 27874 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 26041 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 21521 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 24402 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 19968 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 21270 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 31354 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 38388 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:238: RuntimeWarning: Glyph 35889 missing from current font.\n",
      "  font.set_text(s, 0.0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 26469 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 27874 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 26041 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 21521 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 24402 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 19968 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 21270 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 31354 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 38388 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n",
      "/home/gl/venv0825/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py:201: RuntimeWarning: Glyph 35889 missing from current font.\n",
      "  font.set_text(s, 0, flags=flags)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAR9ElEQVR4nO3de7BdZX3G8e9jIghVEQggEtKg4Ci1VcddtN6qcjE6aqylBW9N64WxU2rV6bQ4drygjlCt2k7RNhVbdKroeM2M2oAB7RXLAW0Vb6R4IYgSTKRlEBH49Y+9oieHfZKdl5y99sn5fmbO7PW+6032L+8k82St9a61UlVIkrSn7tF3AZKkxckAkSQ1MUAkSU0MEElSEwNEktRked8FTNKKFStq9erVfZchSYvGihUr2Lhx48aqWjN335IKkNWrVzMzM9N3GZK0qCRZMarfU1iSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKaGCCSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKaGCCSpCa9BkiSNUm+kWRzkrNG7N8/yYe6/V9IsnrO/lVJbk7yxxMrWpIE9BggSZYB5wFPA44Hnpvk+DnDXgxsr6pjgXcA587Z/3bgMwtdqyTprvo8AjkB2FxV11TVbcCFwNo5Y9YCF3TbHwFOTBKAJM8GvgVcNZlyJUmz9RkgRwHXzmpv6fpGjqmq24GbgEOT3Bv4U+ANu/uSJGckmUkys3Xr1r1SuCRp8V5Efz3wjqq6eXcDq2p9VQ2qanDYYYctfGWStEQs7/G7rwOOntVe2fWNGrMlyXLgIOCHwKOBU5P8OXA/4M4kt1bVXy941ZIkoN8AuRw4LskxDIPidOB5c8ZsANYB/wGcClxSVQU8YceAJK8HbjY8JGmyeguQqro9yZnARmAZ8N6quirJ2cBMVW0Azgfen2QzsI1hyEiSpkCG/6FfGgaDQc3MzPRdhiQtKkmuqKrB3P7FehFdktQzA0SS1MQAkSQ1MUAkSU0MEElSEwNEktTEAJEkNTFAJElNDBBJUhMDRJLUxACRJDUxQCRJTQwQSVITA0SS1MQAkSQ1MUAkSU0MEElSEwNEktTEAJEkNTFAJElNDBBJUhMDRJLUxACRJDUxQCRJTQwQSVITA0SS1MQAkSQ1MUAkSU0MEElSEwNEktSk1wBJsibJN5JsTnLWiP37J/lQt/8LSVZ3/ScnuSLJl7vPp0y8eEla4noLkCTLgPOApwHHA89NcvycYS8GtlfVscA7gHO7/huBZ1bVLwPrgPdPpmpJ0g59HoGcAGyuqmuq6jbgQmDtnDFrgQu67Y8AJyZJVX2xqr7X9V8FHJBk/4lULUkC+g2Qo4BrZ7W3dH0jx1TV7cBNwKFzxvwmcGVV/WSB6pQkjbC87wLujiS/xPC01im7GHMGcAbAqlWrJlSZJO37+jwCuQ44elZ7Zdc3ckyS5cBBwA+79krg48DvVNX/zPclVbW+qgZVNTjssMP2YvmStLT1GSCXA8clOSbJfsDpwIY5YzYwvEgOcCpwSVVVkvsBnwLOqqp/m1TBkqSf6y1AumsaZwIbga8BH66qq5KcneRZ3bDzgUOTbAZeBexY6nsmcCzw2iRf6n4On/AfQZKWtFRV3zVMzGAwqJmZmb7LkKRFJckVVTWY2++d6JKkJgaIJKmJASJJamKASJKaGCCSpCYGiCSpyR4HSJJf6J6kK0lawnYbIEnukeR5ST6V5Abg68D1Sb6a5K1Jjl34MiVJ02acI5BLgQcBrwbuX1VHV9XhwOOBy4Bzk7xgAWuUJE2hcZ7Ge1JV/XRuZ1VtAz4KfDTJPfd6ZZKkqbbbAJkdHkkOBh4A/Bj4dlXdOXeMJGlp2G2AJDkI+APgucB+wFbgXsARSS4D3lVVly5olZKkqTPOKayPAO8DnlBVP5q9I8mjgBcmeWBVnb8A9UmSptQ4p7BO3sW+K4Ar9mpFkqRFYaz7QJIsT5Ju++gkpyZ55MKWJkmaZuPcB/JS4AbgO932JoZvB7wwyZ8ucH2SpCk1zjWQVzC8D+Q+DN8c+ItVdWOSAxm+lvbchStPkjStxgmQ26pqO7A9yeaquhGgqm5JctvClidJmlbjBMgB3fWOewD7ddvpfu61kMVJkqbXOAHyfeDtI7Z3tCVJS9A4y3ifNIE6JEmLzDh3oj9nV/ur6mN7rxxJ0mIxzimsZ3afhwOPBS7p2k8G/h0wQCRpCRrnFNbvASS5CDi+qq7v2kcC/7Cg1UmSptaevJHw6B3h0fkBsGov1yNJWiTGOYW1w6YkG4EPdu3TgM/u/ZIkSYvB2AFSVWcm+Q3giV3X+qr6+MKUJUmaduOswkpVFUAXGHcJjdljJElLw1jvRE/yh0l2ut6RZL8kT0lyAbBuYcqTJE2rcU5hrQFeBHwwyQOB7cABDMPnIuCdVfXFhStRkjSNxlnGeyvwLuBdSe4JrAB+PPfthJKkpWWsi+hJXjuib3bzhqr6mz398iRrgL8ElgHvqapz5uzfn+HrdB8F/BA4raq+3e17NfBi4A7g5VW1cU+/X5LUbtxVWI8BTmf4BN5RLgD2KECSLAPOA04GtgCXJ9lQVV+dNezFwPaqOjbJ6QzfPXJakuO7en4JeADw2SQPrqo79qQGSVK7cW8kvKOq/reqbhr1A7SswDoB2FxV11TVbcCFwNo5Y9YyDCeAjwAndq/WXQtcWFU/qapvAZu730+SNCHjBsjuAqIlQI4Crp3V3tL1jRxTVbcDNwGHjvlrAUhyRpKZJDNbt25tKFOSNMq4p7DumeS+8+wLw2sYU6mq1gPrAQaDgfeqSNJeMm6AXMbw3ejz+UzDd18HHD2rvbLrGzVmS5LlwEEML6aP82slSQtoTx6mmF38tLgcOC7JMUn2Y3hRfMOcMRv4+U2KpwKXdHe8bwBOT7J/kmOA44D/bKxDktRg3COQR7OXV2FV1e1JzgQ2MjwF9t6quirJ2cBMVW0Azgfen2QzsK2rgW7ch4GvArcDf+AKLEmarHED5I6q+t/5diZpurZQVZ8GPj2n77Wztm8FfmueX/tm4M0t3ytJuvv6XIUlSVrE9vlVWJKkhdHnKixJ0iK2J28kbF1tJUnaB/W2CkuStLj1ugpLkrR4uQpLktTEVViSpCauwpIkNXEVliSpiauwJElNXIUlSWriKixJUhNXYUmSmrgKS5LUxFVYkqQmrsKSJDVxFZYkqYmrsCRJTVyFJUlqsjdWYQVXYUnSkuNFdElSEy+iS5KaeBFdktTEi+iSpCZ7ehF9vmsg/7RXqpEkLRpjBUhVvWGhC5EkLS7jXgORJGknBogkqYkBIklq0kuAJDkkycVJru4+D55n3LpuzNVJ1nV9Byb5VJKvJ7kqyTmTrV6SBP0dgZwFbKqq44BNXXsnSQ4BXsfwLvgTgNfNCpq3VdVDgEcCj0vytMmULUnaoa8AWcvw8Sd0n88eMeapwMVVta2qtgMXA2uq6paquhSgqm4DrgRWLnzJkqTZ+gqQI6rq+m77+8ARI8YcBVw7q72l6/uZJPcDnsnwKEaSNEF78krbPZLks8D9R+x6zexGVVXLs7SSLAc+CPxVVV2zi3FnAGcArFq1ak+/RpI0jwULkKo6ab59SX6Q5Miquj7JkcANI4ZdBzxpVnsl8LlZ7fXA1VX1zt3Usb4by2Aw8JldkrSX9HUKawOwrtteB3xyxJiNwClJDu4unp/S9ZHkTcBBzP+OEknSAusrQM4BTk5yNXBS1ybJIMl7AKpqG/BG4PLu5+yq2pZkJcPTYMcDVyb5UpKX9PGHkKSlLFVL56zOYDComZmZvsuQpEUlyRVVNZjb753okqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKaGCCSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKaGCCSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJatJLgCQ5JMnFSa7uPg+eZ9y6bszVSdaN2L8hyVcWvmJJ0lx9HYGcBWyqquOATV17J0kOAV4HPBo4AXjd7KBJ8hzg5smUK0maq68AWQtc0G1fADx7xJinAhdX1baq2g5cDKwBSHJv4FXAmxa+VEnSKH0FyBFVdX23/X3giBFjjgKundXe0vUBvBH4C+CW3X1RkjOSzCSZ2bp1690oWZI02/KF+o2TfBa4/4hdr5ndqKpKUnvw+z4CeFBVvTLJ6t2Nr6r1wHqAwWAw9vdIknZtwQKkqk6ab1+SHyQ5sqquT3IkcMOIYdcBT5rVXgl8Dvg1YJDk2wzrPzzJ56rqSUiSJqavU1gbgB2rqtYBnxwxZiNwSpKDu4vnpwAbq+rdVfWAqloNPB74puEhSZPXV4CcA5yc5GrgpK5NkkGS9wBU1TaG1zou737O7vokSVMgVUvnssBgMKiZmZm+y5CkRSXJFVU1mNvvneiSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKaGCCSpCYGiCSpiQEiSWpigEiSmhggkqQmBogkqYkBIklqYoBIkpoYIJKkJgaIJKmJASJJamKASJKapKr6rmFikmwFvtN3HT1bAdzYdxFTxPnYmfNxV0t9Tm4EqKo1c3csqQARJJmpqkHfdUwL52NnzsddOSfz8xSWJKmJASJJamKALD3r+y5gyjgfO3M+7so5mYfXQCRJTTwCkSQ1MUAkSU0MkH1YklcmuSrJV5J8MMm9MvTmJN9M8rUkL++7zkmZZz5OTHJlki8l+dckx/Zd5yQl+aNuPq5K8oqu75AkFye5uvs8uOcyJ2ae+Xhrkq8n+e8kH09yv36rnB4GyD4qyVHAy4FBVT0MWAacDvwucDTwkKp6KHBhb0VO0C7m493A86vqEcAHgD/rrcgJS/Iw4KXACcDDgWd0AXoWsKmqjgM2de193i7m42LgYVX1K8A3gVf3V+V0MUD2bcuBA5IsBw4Evgf8PnB2Vd0JUFU39FjfpI2ajwLu2+0/qOtbKh4KfKGqbqmq24HPA88B1gIXdGMuAJ7dT3kTN3I+quqirg1wGbCytwqnjAGyj6qq64C3Ad8FrgduqqqLgAcBpyWZSfKZJMf1Week7GI+XgJ8OskW4IXAOf1VOXFfAZ6Q5NAkBwJPZ3h0ekRVXd+N+T5wRF8FTth88zHbi4DPTLyyKWWA7KO689ZrgWOABwC/kOQFwP7Ard2jGf4OeG9/VU7OLubjlcDTq2ol8PfA2/urcrKq6mvAucBFwD8BXwLumDOmGB6l7fN2Nx9JXgPcDvxjH/VNIwNk33US8K2q2lpVPwU+BjwW2NJtA3wc+JWe6pu0UfPxOODhVfWFbsyHGM7RklFV51fVo6rqicB2huf4f5DkSIDuc8mc5pxnPkjyu8AzGF4vWxKBOg4DZN/1XeAxSQ5MEuBE4GvAJ4And2N+ne4fyBIwaj6+ChyU5MHdmJMZztGSkeTw7nMVw+sfHwA2AOu6IeuAT/ZT3eSNmo8ka4A/AZ5VVbf0Wd+08U70fViSNwCnMTzs/iLD8/0HMDwEXwXcDLysqv6rtyInaJ75eDpwNnAnw/9xvqiqrumtyAlL8i/AocBPgVdV1aYkhwIfZvh35DvAb1fVth7LnJh55mMzw1O/P+yGXVZVL+urxmligEiSmngKS5LUxACRJDUxQCRJTQwQSVITA0SS1MQAkaZYkgOSfD7Jsq79yu7pwad17f2S/HP3fC9povxLJ90NSV4PPIbhvSUw/Dd12Tx97El/Vb2e4bOXPlZVdyS5N/CrDJ8W+1HgQ1V1W5JNDO9v8REbmigDRLr7Tq+qHwF074p4xTx9843dVf/zged12+k+59689QngLRggmjBPYUlTKsl+wAOr6tsAVfV/wJeBGYbP7drhKwyPTKSJ8ghEml4rgB/N7qiqtzA82pjdd0eS25LcpwsZaSI8ApGm14+Be405dn/g1gWsRboLA0SaUlW1HViWZJch0j388MbuMfXSxBgg0nS7CHj8bsY8GfjUBGqRdmKASNPtPH7+bo75PA/42wnUIu3Ei+jS3XMD8L4kd3btezB8HeqoPva0v6quTHJpkmVVtdPrZuFnK7U+UVVL5cVgmiK+D0SS1MRTWJKkJgaIJKmJASJJamKASJKaGCCSpCb/D4PkNz+y3BfIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"helloworld\")\n",
    "from main import mydoaconfig\n",
    "import scipy.fft\n",
    "import numpy as np\n",
    "from scipy.linalg import inv, eig, svd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.signal import hilbert, chirp\n",
    "import scipy.constants as constv\n",
    "\n",
    "class Doaest:\n",
    "\n",
    "    def __init__(self, parameter):\n",
    "        self.parameter = parameter\n",
    "        \n",
    "        A = np.random.rand(181, 91, 6) + 1j * np.random.rand(181, 91, 6)\n",
    "        x=[]\n",
    "        y=[]\n",
    "        if parameter.ARRAYTYPE=='ULA':      \n",
    "            x=np.array(parameter.SPACING*0.01*np.linspace(0, parameter.ELEMENTCNT-1, num=parameter.ELEMENTCNT, endpoint=True))\n",
    "            y=np.array([0.,0.,0.,0.,0.,0.])\n",
    "        else:\n",
    "            elementagnle=np.linspace(0, 360, num=6, endpoint=False)/180*constv.pi\n",
    "            x=np.array(parameter.SPACING*0.01*np.cos(elementagnle))\n",
    "            y=np.array(parameter.SPACING*0.01*np.sin(elementagnle))\n",
    "            \n",
    "        pi=constv.pi\n",
    "        l = constv.c/parameter.FREQUENCE/1e9\n",
    "        philist = np.linspace(0, 90, num=91, endpoint=True)\n",
    "        thetalist = np.linspace(0, 180, num=181, endpoint=True)\n",
    "\n",
    "        ii=0\n",
    "        jj=0\n",
    "        for theta in thetalist:\n",
    "            for phi in philist:\n",
    "                tao = x*(((np.sin(phi*pi/180))*(np.cos(theta*pi/180))))+y*(((np.sin(phi*pi/180))*(np.sin(theta*pi/180))));\n",
    "                A[ii][jj] = np.exp(-1j*2*pi*tao/l);\n",
    "                jj=jj+1\n",
    "            ii=ii+1\n",
    "            jj=0\n",
    "        self.A181=A\n",
    "        print(self.A181)\n",
    "        \n",
    "    \n",
    "    def __str__(self):\n",
    "        return (self.parameter.__str__())\n",
    "    \n",
    "    def doaest(self,X):\n",
    "        print(\"doaest\")\n",
    "        R=np.dot(X,(X.conj().T))\n",
    "        a=1/np.trace(R)\n",
    "        R_power = np.linalg.matrix_power(a*R, self.parameter.POWER)\n",
    "        Rbeinv = R_power + np.diag(np.ones(R.shape[0]))\n",
    "        U = []\n",
    "        # Inversion of matrix A\n",
    "        if np.allclose(np.linalg.det(Rbeinv), 0, atol=1e-8):  # Check if A is invertible\n",
    "            print(\"\\nMatrix A is not invertible.\")\n",
    "        else:\n",
    "            U = inv(Rbeinv)\n",
    "            print(\"\\nInversion of A:\\n\", U)\n",
    "        if np.array(U).shape[0]==0:\n",
    "            return -1\n",
    "        print(\"U\\n\",U)\n",
    "        \n",
    "        \n",
    "        PU= np.random.rand(181, 91) + 1j * np.random.rand(181, 91)\n",
    "        for ii in range(181):\n",
    "            for jj in range(91):\n",
    "                #print(self.A181[ii][jj])\n",
    "                #print(U[:,0])                \n",
    "                PU[ii][jj]=(np.dot(self.A181[ii][jj],U[:,0]))\n",
    "                jj=jj+1\n",
    "            ii=ii+1\n",
    "            jj=0\n",
    "        print('pu.shape',PU.shape)\n",
    "        \n",
    "        plt.plot(PU[:,10])\n",
    "        \n",
    "\n",
    "        return 0\n",
    "\n",
    "\n",
    "# Example usage:\n",
    "doaconfig = mydoaconfig(\"./config.ini\")\n",
    "doaestclass = Doaest(doaconfig)\n",
    "\n",
    "print(doaestclass) \n",
    "\n",
    "X=np.array(np.diag([1,1,1,1,1,1]))\n",
    "\n",
    "\n",
    "\n",
    "# array signal parameters\n",
    "fc = 77e9     # carrier frequency\n",
    "c = 3e8      # light speed\n",
    "M = 6        # array sensor number\n",
    "N = 400       # snapshot number\n",
    "wavelength = c / fc  # signal wavelength\n",
    "d = 0.5 * wavelength  # inter-sensor distance\n",
    "SNR = 10\n",
    "DOA = np.array([0,10])\n",
    "\n",
    "                    \n",
    "array_geom = np.expand_dims(np.array(np.arange(M)), axis=-1) * d\n",
    "\n",
    "signal_i = 10 ** (SNR / 20) * (np.random.randn(K, N) + 1j * np.random.randn(K, N))\n",
    "noise = np.random.randn(M, N) + 1j * np.random.randn(M, N)\n",
    "\n",
    "phase_shift = 2 * np.pi * array_geom / wavelength * np.sin(DOA / 180 * np.pi)\n",
    "a_i = np.exp(-1j*phase_shift)\n",
    "array_signal = np.matmul(a_i, signal_i) \n",
    "print(\"array_signal_shape:\",array_signal.shape)\n",
    "    \n",
    "doaestclass.doaest(array_signal)\n",
    "\n",
    "\n",
    "K = len(DOA)\n",
    "w,V = np.linalg.eig(np.dot(array_signal,array_signal.conj().T))\n",
    "w_index_order = np.argsort(w)\n",
    "V_noise = V[:,w_index_order[0:-K]]\n",
    "noise_subspace = (np.matmul(V_noise, np.matrix.getH(V_noise)))\n",
    "doa_search = np.arange(-90,90)\n",
    "p_music = np.zeros(len(doa_search))\n",
    "K = len(DOA.reshape(-1))\n",
    "array_geom = np.expand_dims(np.array(np.arange(M)), axis=-1) * d\n",
    "phase_shift = 2 * np.pi * array_geom / wavelength * np.sin(DOA / 180 * np.pi)\n",
    "a = np.exp(-1j*phase_shift)\n",
    "\n",
    "for doa_index in range(len(doa_search)):\n",
    "\n",
    "    K = len(DOA)\n",
    "    w,V = np.linalg.eig(np.dot(array_signal,array_signal.conj().T))\n",
    "    w_index_order = np.argsort(w)\n",
    "    V_noise = V[:,w_index_order[0:-K]]\n",
    "    noise_subspace = (np.matmul(V_noise, np.matrix.getH(V_noise)))\n",
    "    doa_search = np.arange(-90,90)\n",
    "    p_music = np.zeros(len(doa_search))\n",
    "    K = len(DOA.reshape(-1))\n",
    "    array_geom = np.expand_dims(np.array(np.arange(M)), axis=-1) * d\n",
    "    phase_shift = 2 * np.pi * array_geom / wavelength * np.sin(doa_search[doa_index] / 180 * np.pi)\n",
    "    a = np.exp(-1j*phase_shift)\n",
    "\n",
    "\n",
    "    p_music[doa_index]=1/np.matmul(np.matmul(a.T.conjugate(),noise_subspace),a).reshape(-1)[0]\n",
    "    \n",
    "    \n",
    "p_music = p_music/np.max(p_music)\n",
    "p_music = 10*np.log10(p_music)\n",
    "plt.figure()\n",
    "plt.plot(doa_search,p_music)\n",
    "plt.ylabel(\"归一化空间谱(dB))\", fontproperties=\"SimSun\")\n",
    "plt.xlabel(\"来波方向(°)\", fontproperties=\"SimSun\")\n",
    "plt.show()\n",
    "\n",
    "                    \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv0825",
   "language": "python",
   "name": "venv0825"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
